DOE/ER/53223-34 
MF-115 


Courant  Institute  of 
Mathematical  Sciences 

Magneto-Fluid  Dynamics  Division 


Magnetic  Surfaces  and  Neoclassical 
Transport  in  Stellarators 


Kam-Chuen  Ng 


S  2  S 

cn  3   0)  lO  P? 

n-  rD  iX)  -  G 
0)   O   3 

M  M  n-  cj 

S  S  of  ?  Supported  by  U.S.  Department  of  Energy 
g  H.m  I  '    Grant  No.  DE-FG02-86ER53223;  and  NSF 
m  ^  ;;  S  M  Grant  No.  DMS-8320430. 

rr  O        ui 

S  m  Plasma  Physics 

w  o'  June  1987 

O    Dj 

1-1  M 

rr 

(-•• 

3 


NEW  YORK  UNIVERSITY 


UNCLASSIFIED 


DOE/ER/53223-34 

UC-20 

Plasma  Physics 


Courant  Institute  of  Mathenatical  Sciences 
New  York  University 


MAGNETIC  SURFACES  AND  NEOCLASSICAL  TRANSPORT 

IN  STELLARATORS 

Kam-Chuen  Ng 

June  1987 


This  work  was  supported  by  the  U.  S.  Department  of  Energy 
Grant  No.  DE-FG02-86SR53223;  and  NSF  Grant  No.  DMS-8320436; 
some   of  the  computations   were  performed  at  the   San  Diego 

Supercomputer  Center. 


UNCLASSIFIED 


DISCLAIMER 


This  report  was  prepared  as  an  account  of  work  sponsored 
by  an  agency  of  the  United  States  Government.   Neither 
the  United  States  Government  nor  any  agency  thereof,  nor 
any  of  their  employees,  makes  any  warranty,  express  or 
implied,  or  assumes  any  legal  liability  or  responsibility 
for  the  accuracy,  completeness,  or  usefulness  of  any 
information,  apparatus,  product,  or  process  disclosed,  or 
represents  that  its  use  would  not  infringe  privately 
owned  rights.   Reference  herein  to  any  specific  commercial 
product,  process,  or  service  by  trade  name,  trademark, 
manufacturer,  or  otherwise,  does  not  necessarily  constitute 
or  imply  its  endorsement,  recommendation,  or  favoring  by 
the  United  States  Government  or  any  agency  thereof.   The 
views  and  opinions  of  authors  expressed  herein  do  not 
necessarily  state  or  reflect  those  of  the  United  States 
Government  or  any  agency  thereof. 


Printed  in  U.S.A. 

Available  from 

National  Technical  Information  Service 
U.S.  Department  of  Commerce 
5285  Port  Royal  Road 
Springfield,  VA  22161 


ABSTRACT 

In   this  thesis  we  shall  study  the  structure  of  a  stellarator  field 

and   the   confinement  of  a  high  temperature  plasma  in  toroidal  geometry. 

A  field  line   tracing  program  is  developed  to  explore  the  structure  of 

magnetic   fields   on   a   fine   scale  so  as  to  explain  anomalous  electron 

transport.    The   model   magnetic   field  we  choose  has  a  simple  analytic 

representation  which   is   easy  to  compute ,  so  that  we  can  integrate  the 

field   lines  to  a  high  accuracy.   On  the  CRAY  computer,  integration  over 

one   field  period   requires   only   SO^s .   In  a  typical  case  most  of  the 

magnetic  surfaces  are  well  behaved  on  the  scale  of  the  gyroradius  of  the 

electron,  p    ,   even  when  the  magnetic  field  has  no  2d  symmetry.   Island 

chains   or   stochastic   regions   are   formed  in  the  vicinity  of  magnetic 

surfaces   with   rational  rotational  transform  t,   -  ~    .      We  show  that  the 

m 

island  width  w  decays  exponentially  with  m.  We  have  calculated  an 
example  where  the  island  width  w  is  less  than  the  gyroradius  of  the 
electron  p  for  m  greater  than  17,  which  means  that  higher  order 
islands  will  not  affect  the  electron  transport. 

Our  results  suggest  that  the  anomalous  electron  transport  observed 
in  experiments  may  be  due  to  the  presence  of  an  ambipolar  electrostatic 
potential  $.  We  investigate  this  hypothesis  by  computing  the  guiding 
center  orbits  of  the  electrons  and  estimating  island  widths  of  the  drift 
surfaces  that  are  swept  out.  We  find  that  with  a  small  electric 
potential  depending  on  the  toroidal  and  poloidal  angles,  the  drift 
surface  island  width  w  is  an  order  of  magnitude  larger  than  that  of  the 
magnetic   surfaces  and  decays  exponentially  at  a  slower  rate.   Since  the 


11 


drift  step   size   is   the   maximum  of  p      and  w,  in  this  case  electron 

e 

transport,  which  scales  like  the  square  of  the  step  size,  is  enhanced. 
Numerical  results  confirm  that  resonance  of  the  electric  potential  *  at 
the  rational  magnetic  surfaces  may  cause  anomalous  electron  transport. 
The  conclusion  is  substantiated  by  a  calculation  of  the  confinement  time 
using  the  Monte  Carlo  method.  Introduction  of  an  electric  potential  as 
small  as  two  percent  of  the  temperature  may  bring  down  the  confinement 
time  by  an  order  of  magnitude. 
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I.   GENERAL  DISCUSSION  OF  THE  PROBLEM 


A  principal  goal  in  controlled  thermonuclear  fusion  research  is  to 
confine  high  temperature  plasma  by  a  strong  magnetic  field  in  toroidal 
geometry  long  enough  for  the  nuclear  reactions  to  produce  more  energy 
than  is  consumed.  The  confinement  of  plasma  in  such  geometry  depends  on 
the  existence  of  a  family  of  nested  magnetic  flux  surfaces.  According 
to  the  ideal  magnetohydrodynamic  (MHD)  model  of  a  perfectly  conducting 
fluid,  plasma  cannot  flow  across  the  magnetic  field.  Therefore  the 
existence  of  nested  flux  surfaces  guarantees  the  confinement  of  the 
plasma.  However,  in  practice  the  situation  is  not  that  simple.  The 
existence  of  nested  flux  surfaces  has  only  been  proved  in  cases  with 
two-dimensional  symmetry.  In  general  the  experimental  devices  have  no 
symmetry  and  we  do  not  have  perfect  confinement.  The  plasma  in  the 
container  will  eventually  diffuse  to  the  wall  due  to  the  effects  of 
resistivity  and  finite  gyroradius . 

In  experiments,  the  transport  of  electrons  is  found  to  be  one  or 
two  orders  of  magnitude  larger  than  that  predicted  by  neoclassical 
theory  [16].  The  existing  analytic  theory  does  not  fully  explain  this 
phenomenon.  Recently  Bauer,  Betancourt  and  Garabedian  [5]  developed  a 
Monte  Carlo  code  to  calculate  neoclassical  plasma  transport.  The  code 
couples  the  result  of  a  three-dimensional  MHD  equilibrium  code  with  a 
particle  transport  calculation.  The  numerical  results  suggest  that 
anomalous  electron  transport  may  be  due  to  small  resonant  terms  in  the 
electric  potential.  It  is  our  goal  to  give  an  independent  study  of  the 
problem  by  estimating  the  island  width  of  the  magnetic  surfaces  and  the 
drift  surfaces   numerically.    We   find  that   introduction  of  a  small 


electric  potential  does  indeed  enhance  the  island  width  significantly. 
This  means  that  anomalous  transport  of  electrons  may  be  due  to  the 
resonance  of  the  electric  potential.  Also  our  numerical  results 
concerning  the  island  width  of  drift  surfaces  imply  that  transport 
associated  with  these  islands  is  independent  of  the  mass  of  the 
particle.  This  means  that  electron  transport  and  ion  transport  are  of 
the  same  order  of  magnitude,  which  has  a  bearing  on  charge  neutrality. 
In  the  following  we  shall  discuss  these  important  issues  in  detail. 

1 .   Nested  flux  surface  hypothesis 

If  one  follows  a  field  line  around  the  torus  many  times  and  if  the 
field  line  remains  within  the  torus,  there  are  three  possibilities  for 
the  behavior  of  the  field  line  depending  on  the  degree  of  symmetry  of 
the  magnetic  field.  The  first  possibility  is  that  the  field  line  closes 
itself  after  a  finite  number  of  field  periods.  The  second  possibility 
is  that  the  field  line  covers  a  surface  ergodically.  In  this  case  we 
have  a  magnetic  flux  surface.  The  third  possibility  is  that  the  field 
line  forms  a  stochastic  region  in  space.  The  last  case  is  the  most 
general.  We  know  that  we  have  nested  magnetic  surfaces  everywhere  when 
the  magnetic  field  has  2d  symmetry.  The  existence  of  magnetic  surfaces 
in  general  is  related  to  the  divergence  free  property  of  the  magnetic 
field  [12] . 

Since  the  magnetic  field  is  divergence  -  free ,  the  Poincare  mapping 
defined  by  following  the  field  lines  around  the  torus  once  is  measure 
preserving,  so  the  Kolmogorov-Arnold-Moser  (KAM)  theory  of  Hamiltonian 
systems  can  be  applied.  The  existence  of  magnetic  surfaces  in  certain 
cases  is  guaranteed  by  the  KAM  theory,  which  we  will  describe  briefly  in 


the  context  of  plasma  confinement  in  a  torus.  The  KAM  theory  is 
concerned  with  invariant  curves  of  a  Hamiltonian  system  under 
perturbation.   Let  us  define  these  concepts  in  more  detail. 

Let  T  be  the  transformation  of  a  plane  cross  section  C  of  the 
torus  onto  itself  defined  by  following  the  field  lines  once  around  the 
torus.  If  X  is  a  point  on  C,  let  Tx  stand  for  the  image  point  of  x  on 
C.  T  X  describes  the  successive  iterations  of  the  point  x-  A  curve  T 
on  C  is  invariant  if  for  any  point  on  the  curve  F  its  image  under  the 
transformation  T  is  also  on  that  curve.  We  consider  the  iterates  of  a 
single  point  on  a  closed  invariant  curve  T  and  introduce  the  expression 


1     1  N 
T^  lim  ^  •£.   e 

N-»co   n-1 


where  d  is  the  angular  displacement  of  the  nth  iterate  of  the  mapping  T 
with  respect  to  a  fixed  point  inside  the  closed  invariant  curve  T.  If  t 
exists  then  its  value  is  independent  of  the  choice  of  the  fixed  point 
lying  inside  the  closed  curve  and  of  the  initial  point  of  the  invariant 
curve  r.  The  number  t(r)  is  called  the  rotational  transform  of  the 
invariant  curve  T.  The  property  that  B  is  divergence  free  allows  us  to 
apply  the  KAM  theory  here,  and  more  specifically  Moser's  theorem  on 
measure-preserving  mappings  [21].  Since  B  is  divergence  free,  the 
divergence  theorem  shows  chat  the  transformation  T  is  measure-preserving 
in  the  sense  that  it  leaves  invariant  the  magnetic  flux 
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where  S  is  any  region  on  the  cross  section  C. 


Suppose  T„  is  a  measure -preserving  mapping  which  possesses  a  family 
of  nested  closed  invariant  curves  r  bounded  inside  by  T  and  outside  by 
r,  .   Assume  that  the  associated  rotational  transforms  t(r)  in  the  range 


0  <  t(r^)  <   tin   <   t(rjj) 


have  shear,  i.e.  t'  >  0,  where  the  differentiation  is  with  respect  to 
some  variable  that  defines  the  invariant  curves.  Moser  [21]  proves  that 
for  a  given  £  >  0  there  is  a  6  >  0  such  that  if  a  sufficiently  smooth 
measure -preserving  mapping  T^  satisfies 


MTq  -  t^II  <  « 


and  if  cj  is  any  number  in  the  interval 


t(r  )  +  £  <(j<  t(r,) 

a  ^  b 


which  cannot  be  closely  approximated  by  rationals  p/q  in  the  sense  that 


q    3 
q 


then  T,  will  have  a  closed  invariant  curve  F  with  rotational  transform 
8  u> 

u).  Here  the  positive  number  6  depends  on  t  and  also  on  the  smoothness 
of  the  measure -preserving  mapping  T^ .  In  fact  there  exist  infinitely 
many  closed  invariant  curves,  since  the  set  of  the  numbers  cj  not  closely 
approximated  by  the  rationals  has  positive  measure  close  to  the  measure 
of  the  interval  [t(r  ),t(r,)]  if  «  is  small. 

a,  D 


The  above  theorem  says  that  although  we  cannot  recover  all  the  flux 
surfaces,  after  performing  the  perturbation  T,  we  still  have  infinitely 
many  of  them.  However,  when  the  amplitude  of  the  perturbation  increases 
the  theory  breaks  down.   We  shall  investigate  this  case  numerically. 

An  analysis  of  the  orbit  equations  of  a  charged  particle  shows  that 
magnetic  field  lines  are  the  lowest  order  approximation  for  the  orbits 
of  the  guiding  center.  Therefore  if  the  gyroradius  of  the  charged 
particle  is  small  and  the  magnetic  flux  surfaces  exist,  we  can  apply  the 
KAM  theory  to  conclude  that  nearby  drift  surfaces  exist  and  hence  good 
confinement  of  the  charged  particles  is  expected. 

In  the  equilibrium  and  stability  code  BETA,  which  is  based  on  the 
variational  principle  of  ideal  MHD,  it  is  assumed  that  a  family  of 
nested  magnetic  surfaces  exists  [2].  To  get  a  qualitative  assessment 
concerning  the  adequacy  of  the  equilibrium  code,  Marcal  [18]  developed 
a  method  of  line  tracing  that  was  based  on  a  divergence- free 
approximation  to  the  numerical  equilibrium  generated  by  the  code.  Due 
to  the  crudeness  of  the  mesh,  small  islands  comparable  in  size  to  the 
gyroradius  of  the  ions  or  electrons  cannot  be  detected  using  his  code. 
In  this  paper  we  introduce  a  special  analytical  representation  of  the 
magnetic  field  and  we  are  able  in  the  simplest  cases  to  detect  small 
islands  on  the  length  scale  of  the  gyroradius  p      of  the  electron. 

The  relationship  of  the  hypothesis  of  nested  surfaces  to  the  KAM 
theory  of  Hamiltonian  systems  has  been  discussed  by  Garabedian  [6,10]. 
Even  though  there  may  be  stochastic  regions  where  no  flux  surfaces  are 
swept  out  by  the  field  lines,  the  KAM  theory  predicts  that,  under 
favorable   conditions,   an   infinite   number   of  nested   toroidal   flux 


surfaces  will  still  remain  and  will  predominate  over  the  smaller 
stochastic  regions.  The  finite  number  of  nested  flux  surfaces  implicit 
in  the  computational  grid  used  by  the  equilibrium  code  may  be  thought  of 
as  constituting  a  selection  from  this  infinite  set  of  nested  surfaces. 
However,  it  is  necessary  to  have  independent  ways  of  testing  for  the 
break-up  of  surfaces.  Our  simple  analytical  model  enables  us  to  detect 
small  islands  and  to  study  the  regularity  of  the  magnetic  surfaces  on 
the  length  scale  of  the  gyroradius  of  the  electron.  It  is  found  that 
most  of  the  surfaces  in  a  typical  Heliac  case  behave  regularly  on  this 
fine  scale  (cf.  Chapter  3).  Therefore  the  assumption  of  nested  surfaces 
is  not  unreasonable  if  the  grid  size  is  not  too  small. 

2 .   Transport  and  the  electric  field 

When  plasma  has  a  density  gradient,  it  tends  to  diffuse  towards  the 
regions  of  low  density.  The  central  problem  in  controlled  fusion 
research  is  to  impede  the  rate  of  diffusion  by  means  of  a  magnetic 
field.  The  classical  theory  of  diffusion  of  a  plasma  assumes  that  the 
collisions  between  ions  and  electrons  are  so  rapid  that  a  resistive 
magnetohydrodynamic  model  can  be  used  to  describe  the  behavior  of  the 
plasma.  However,  if  the  density  of  the  plasma  is  very  low  collisions 
are  not  frequent.  A  fluid  model  of  the  plasma  is  not  adequate  to 
explain  its  behavior  and  the  kinetic  theory  becomes  a  more  suitable  tool 
to  explain  the  phenomena  observed  in  experiments.  In  the  regime  of  low 
collisionality ,  the  mean  free  path  of  the  charged  particles  is  much 
greater  than  the  characteristic  length  of  the  device.  The  local 
properties  of  the  magnetic  field  involving  its  gradient  and  curvature 
cause  the  particles  to  drift  across  the  field  lines.  This  effect  may 
increase   the   diffusion   significantly.    In   the   classical   theory,  a 


decrease  in  collisionality  results  in  a  decrease  of  the  diffusion  rate. 
However,  in  actual  cases,  the  diffusion  is  enhanced  by  the  low 
collisionality.  In  particular,  some  of  the  particles  become  trapped  in 
magnetic  wells  produced  by  the  external  helical  windings  in 
stellarators .  These  particles  have  a  relatively  small  component  of 
velocity  parallel  to  the  direction  of  the  magnetic  field,  and  they  drift 
across  the  magnetic  field  lines  as  they  bounce  back  and  forth.  This  is 
a  consequence  of  the  lack  of  symmetry  in  a  stellarator.  To  simulate 
these  effects  one  uses  a  drift  kinetic  equation  with  a  collision 
operator  to  evaluate  the  diffusion.  In  our  work  a  Monte  Carlo  code 
based  on  the  drift  kinetic  equation  has  been  written  to  evaluate  the 
transport  of  the  plasma  (cf.  Chapter  V). 


For   equal   ion  and  electron  temperatures  the  collision  time  scales 

ste^ 
,V2 


1/2 
like   M    ,   the   square   root  of  the  mass  of  the  particle,  and  the  step 


size   is  proportional   to   the   gyroradius ,  which  is  of  the  order  M 

Therefore   the   transport,  which  scales  like  the  square  of  the  step  size 

1/2 
over   the   time   step,   is   proportional   to   M    .   This  means  that  the 

electrons   are   better   confined  by  a  corresponding  factor.   However,  in 

experiments   the   electron   transport   is   found  to  be  one  or  two  orders 

higher  than  that  predicted  by  the  neoclassical  theory  when  the  effect  of 

the   electric   potential   is   not  included.    With  the  introduction  of  a 

small   electric  potential,   islands   form.    The  effect  of  a  large 

collection   of   these   islands  in  the  drift  surfaces  makes  the  step  size 

scale   like   M  '         (cf.  Section  II. 4).   This  brings  down  the  confinement 

time  of  the  electrons. 


The   conventional  form  of  the  electric  potential  does  not  depend  on 
the   toroidal   and  poloidal   angles.   In  our  model,  which  is  consistent 


with  a  generalized  Ohm's  law  proposed  by  Boozer,  the  electric  potential 
does  depend  on  the  toroidal  and  poloidal  angles  (cf.  Chapter  IV).  This 
form  of  the  electric  field  is  believed  to  be  the  main  factor  that  causes 
resonance  and  enhances  the  electron  transport.  In  Chapter  V  we  present 
a  Monte  Carlo  calculation  of  neoclassical  transport  based  on  our  simple 
analytical  representation  of  the  magnetic  field.  The  results  do  suggest 
that  anomalous  transport  of  electrons  may  be  due  to  small  resonant  terms 
in  the  electric  potential. 
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II.  TWO-DIMENSIONAL  EQUILIBRIUM 


In  this  chapter  we  shall  give  a  brief  description  of  a  variational 
principle  which  characterizes  the  equilibrium  solution  of  the  partial 
differential  equations  of  ideal  magnetohydrodynamics .  Based  on  this 
variational  principle,  a  numerical  procedure  to  study  the  equilibrium 
and  the  stability  of  a  plasma  without  two-dimensional  symmetry  has  been 
developed  by  Bauer,  Betancourt  and  Garabedian  [2,3].  Their  algorithm 
makes  use  of  an  accelerated  method  of  steepest  descent  associated  with 
the  variational  principle.  We  specialize  the  problem  of  finding  an 
equilibrium  solution  to  the  case  of  two-dimensional  helical  symmetry. 
The  solution  of  the  two-dimensional  problem  can  serve  as  an 
initialization  to  a  more  complicated  three-dimensional  problem  so  as  to 
speed  up  the  calculation.  Moreover,  when  the  problem  has  2d  symmetry 
the  2d  code  gives  a  more  efficient  and  accurate  result.  The  results 
from  the  2d  helical  code  are  compared  with  exact  solutions  and  with 
those  of  the  general  3d  code.  They  have  a  bearing  on  our  investigation 
of  magnetic  surfaces  in  two  and  three  dimensions. 

1 .   Variational  principle 

Toroidal   equilibrium   of   a  plasma   can  be  described  by  the  ideal 
magnetostatic  equations 


(2.1)  J  X  B  -  V  p 

V  •  B  -  0 


where  J  -  V  x  B  is  the  current  density,  B  is  the  magnetic  field  and  p  is 
the  scalar  pressure.  The  system  of  partial  differential  equations  (2.1) 
is  of  a  nonstandard  type,  for  it  has  both  real  and  complex 
characteristics.  The  equilibrium  and  stability  of  the  solution  of  (2.1) 
can  be  studied  by  a  classical  variational  principle  for  the  potential 
energy  [15] 

E  -  //J  [  2  ^^  +  :^  1  dxdydz 

subject  to  flux  constraints  and  the  condition  V  •  B  -  0.  The  fluid 
pressure  p  is  related  to  the  density  p  by  p  -  p^ ,  where  7  is  the  ratio 
of  specific  heats. 

We   integrate   V  •  B  -  0  by   representing  B  as  a  cross  product  of 
gradients 
(2.2)  B  -  Vs  X  V0 

where  s  and  0  are  flux  functions.  Here  s  is  assumed  to  be  single -valued 
and  s  =  constant  is  supposed  to  be  a  nested  family  of  toroidal  flux 
surfaces,  so  that  we  can  use  s  as  an  independent  variable.  The  values 
of  s  range  from  0  to  1;  s  =-  0  is  the  magnetic  axis  and  s  =•  1  is  the 
boundary  surface.  On  the  flux  surfaces,  let  u  and  v  be  angle  variables 
with  unit  periods  in  the  poloidal  and  toroidal  directions,  respectively. 
The  flux  function  i/>  has  periods  in  the  poloidal  and  toroidal  directions 
whose  quotient  is  the  rotational  transform  t  and  depends  on  s  only. 
With  appropriate  constraints  on  the  toroidal  and  poloidal  fluxes  and  on 
the  mass  within  each  flux  surface,  the  Euler  equations  of  the 
variational  problem  are  the  magnetostatic  equations  (2.1),  and  the 
potential  energy  E  is  minimized  when  p  and  p  are  functions  of  s  alone 
[2].   This  allows  us  to  impose  an  additional  ergodic  constraint  p  -  p(s) 
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on  the  variational  problem.  The  additional  constraint  and  the 
representation  (2.2)  of  B  reduce  the  system  by  eliminating  two 
equations,  because  B  is  automatically  divergence  free  and 

(2.3)  B  •  Vp  -  0 

is  satisfied.  Integration  of  the  relation  (2.3)  avoids  a  difficulty 
with  the  real  characteristics  of  the  system,  which  prevent  us  from 
formulating  a  well  posed  boundary  value  problem  for  a  system  of  elliptic 
partial  differential  equations.  This  reduction  also  allows  us  to  choose 
z  -  2»rv,  thus  diminishing  the  number  of  unknown  functions  by  one. 

The  crux  of  the  numerical  computation  of  magnetohydrodynamic 
equilibria  is  to  formulate  a  discrete  variational  principle  for  the 
potential  energy  E.  This  is  accomplished  by  writing  down  a  second-order 
accurate  approximation  of  E  based  on  a  rectangular  mesh  over  the  cube  0 
<  s  <  1,  0  <  u  <  1,  0  <v<l.  Setting  the  derivatives  of  the 
approximation  of  E  with  respect  to  the  nodal  values  of  the  unknowns 
equal  to  zero,  we  obtain  finite  difference  equations.  An  iterative 
scheme  to  solve  the  finite  difference  equations  is  then  developed  by 
formulating  a  suitable  finite  difference  approximation  to  the 
artificially  time -dependent  partial  dif  fei.ential  equations  for  paths  of 
steepest  decent. 

2 .   Helical  symmetry 

We  choose  s,  u  and  v  as  independent  variables  to  minimize  the 
potential  energy  E.  All  physical  quantities  are  computed  in  the  s,  u,  v 
coordinate   system.    The   choice  of  these  coordinates  has  the  advantage 
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that  the  constraints  and  the  boundary  conditions  are  easily  incorporated 
into  the  variational  principle,  which  is  important  for  sound  numerical 
analysis  and  computing. 

Now  consider   the   straight  helically   symmetric  case.   The  point 
(x.y.z)  in  rectangular  coordinates  is  represented  by  the  formulas 


x+iy  -  {XQ(v)  +  iyQ(v)+R(s.u,v)[Xj^(u)+iyj^(u)-XQ(v)-iyQ(v)])e^^''*^ 
z  -  2irv 


where  i  -  /T  and  k  is  an  integer.  Here  x+iy  -  [xQ(v)+iyQ(v)  Je^^"^  is 
the  equation  for  the  magnetic  axis  ,  x+iy  -  [x,  (u)+iy^(u)  le^^*^*^  is  the 
boundary  wall,  and  the  radial  variable  R  is  0  when  s  -  0  and  1  when 
s  -  1.  Since  s  is  one  of  the  independent  variable  and  z  -  27rv,  the  x- 
component  of  the  magnetic  field  B  is  given  by 


g  _  a(s.V».x)  _  d(^.x)    1 
x   a(x,y,z)     a(u,v)  D 


where  the  Jacobian  D  of  the  transformation  is 


Q  _  a(x.y,z)  _  2^  3(x.y) 
a(s,u,v)       5(s,u) 


Similarly  we  have 


„     3(^,y)  1     „     d(i>,z)   1 
y    a(u,v)  D  '     z    a(u,v)  D 


Therefore 
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B     -  -^   {V>   (x  +y  +z   )    -    2^-  V"   (x  X  +y  y  +z   z    )   +  V   (x  +y  +z    )  ) 
2        u     V  ''v     v'  '^u'^v^    u  V  -'u-'v     u  v'         w^   u  ■'u     u 


Since  x,   and  y^   are   functions  of  u  only,   D  does  not  depend  on  v 
explicitly. 


Now  let  X+iY  =  XQ+iyQ+R[x^+iy^-XQ-iyQ]  .   Then  x+iy  -  (X+iY)e^^''^ 
and  we  have 


x^  +  y^  -  (27rk)^(X^  +  Y^)  +  4»rk(XY  -  X  Y)  +  X^  +  Y^ 

V    ■'v  V    V  V     V 

XX  +  y  y  =  2;rk(XY  -  YX  )  +  (X  X  +  Y  Y  ) 

uv^u^v  u  u      uv    uv 

2  2  2    2 

x""  +  y'^  -  X'^  +  Y'^ 

u  ^u  u    u 


From  the  definition  of  X  +  iY,  v  does  not  appear  explicitly,  nor  does  it 

2    2  2    2 

appear  explicitly  in  the  expressions  x  +y   ,xx  +yy  and  x  +  y  . 

t^"^         r         J  r  V-^VUVUV         UU 

In  general  the  flux  function  V  is  given  by 


V-  -  V(s,u,v)  -  -F^(s)u  +  F^(s)v  +  A(s,u,v) 


where   A   is   a   single -valued   function  and  the  quotient  F'/F'   is  the 
rotational  transform  t .   Now  the  potential  energy  E  can  be  written  as 


E  -  ;/J  ^  f^l'^S  -  2Ai2Vv  ^  A22^J)dsdudv  .  -^^   ^(j-j-p^l!)^'^  "' 


where 


A,,    -    (27r)^(l+k^(X^+Y^))    +   47rk(XY    -X  Y)    +   X^    +   Y^ 
W         ^       '     ^  ^  '  '  vv  v  V 


A,  „    =   27rk(XY      -    YX   )    +    (X  X      +  Y  Y   ) 
12  u  u  u  V  u  v' 
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2    2 


Since  v  does  not  appear  explicitly  in  the  integrand,  we  can  seek  a 
solution  il>,  R,  x„,  y^  which  is  independent  of  v  when  minimizing  the 
energy  E.   A  computer  code  has  been  written  to  perform  this  computation. 

3 .   Comparison  with  an  exact  solution 

In  order  to  verify  an  exact  solution  to  our  problem,  we  make  a 
variation  of  the  potential  energy  E  with  respect  to  i>,  R,  Xq  and  y^ . 
The  resulting  Euler  equations  are 


2 


L2(R)    -   Inmj^ij-  +  p) 


7  2 

(27rk)>^X    -    27rkV>^Y^t  ^      2>rkV.^Yt    +  X^i 

+  ('^r^o^'                    5  "  i^^             D              ^' 

(27T-k)^0  Y    -    2jrkV>  X    i  -      27rk0  Xi    +  Y    t 

.,______u u   u  d     ,  u u      ,  „ 

-^  (yi-yo^^ D  -  du^        D         )^-° 


2                                 (2;rk)'^\&^X    -    2jrkV'^Y^t 
^Cxq)    -  JJ    (    27rRR^(y^)^(—  +  p)    +    (1-R)[ ^—5 ^^"^l 

2nk\l>  Yt    +  X   t^ 

-    R    ( ^ ^— ))dsdu  -  0 

u  D 


2                                 (27rk)^V^Y    -    2;rkV.   X^t 
L4(yo)    -  //    (-2'rRR^(x^)^(f-   +   p)    +    (1-R)[ ^ ^"^ 

27rk0   Xi    -    Y    (.^ 

+  R    ( ^-^^ ^^— ))dsdu  -  0 

u  D 
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Here  we  have  put  Fp(s)   =  1  ,  so  that  F'(s)  =  t(s)  is  the  rotational 
transform.  The  exact  solution,  given  by  Schwenn  [23],  is  as  follows.  The 

cross   section  x,+iy, -Xq- iyQ  =  e     is  circular,  k-1,  t/>  -1,  i-2, 

2 
p  -   (l-s)/jr  ,   and  ^q  "  ^0  "  ^/^'   ^®  ^®®  that  this  solution  satisfies 

L,(0)   -  LjCR)  -  L^(Xq)  =  L, (yQ)  -  0,  and  the  2d  computer  code  gives  the 

same   result.    We   have   also   compared  the  results  of  the  2d  code  with 

corresponding  data  from  the  3d  code  [3].   Good  agreement  is  obtained. 


4.   Drift  step  size 

In  this  section  we  employ  the  MHD  model  to  derive  expressions  for 
the  drift  step  size  that  explain  the  breakdown  of  rational  surfaces  with 
the  introduction  of  a  small  electric  potential.  We  start  by  restating 
the  magnetostatic  equations. 

The  partial  differential  equations  for  MHD  equilibrium  are 

(2.5)  V  .  B  -  0 

(2.6)  J  X  B  -  Vp 

(2.7)  J  -  V  X  B 

where  J  is  the  current  density  and  p  is  the  fluid  pressure.  A  basic 
assumption  in  the  formulation  is  that  there  exists  a  nested  set  of  flux 
surfaces  s  —  const,  with  the  pressure  restricted  to  be  a  function  p  - 
p(s)  of  s  alone.   It  is  known  that 

(2.8)  B  -  Vs  X  V0 
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can  be   represented   in   terms  of  two  flux  functions  [1].   The  relation 
(2.6)  implies  that 


Vs  -  0 


To  prove  this,  taking  a  scalar  product  of  equation  (2.6)  with  J  we 
obtain  J  •  Vp  -  0.  If  p'(s)  r*  0  then  the  result  follows  immediately.  If 
p'(s)  -  0  then  we  have  J  x  B  -  0  which  means  that  J  is  a  scalar  multiple 
of  the  vector  B.   Since  B  •  Vs  -  0 ,  we  have  J  •  Vs  -  0 . 

The  relation  J  •  Vs  -  0   implies  that  the  line  integral 


4,{s    ,u,v)  =  J  B'di  -  /  d0 
°        C        C 


is   path   independent,  where  C  is  a  curve  joining  two  points  (u  ,v  )  and 
(u,v)  on  any  surface  s  -  s  .   Here  dZ    is  a  directed  line  element.   Thus 
B  -  V^   is  a  scalar  multiple  of  Vs  and  therefore  B  can  be  represented  as 

(2.9)  B  -  Vai  +  ^Vs 

Let  u  and  v  be  the  angle  variables  defined  in  the  Section  1  that 
have  poloidal  periods  1  and  0  and  toroidal  periods  0  and  1  respectively. 
On  each  flux  surface  s  =  const.  the  functions  0  and  4>  have  a 
representation  [5] 

0  =•  u  -  tv  +  l/).. 

4>   -  I„(s)u  +  I„(s)v  +  4>, 
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whe 


re  0,   and  4),    are   single-valued.   I„   is  the  toroidal  current,  which 


is   zero   in  the  stellarator  case,  and  !_  is   the  poloidal  current.   We 
introduce  renormalized  variables  ^  and  <l>   by  the  affine  transformation 


yl)   -   i>   +    10/1 


where  i>  and  J  are  supposed  to  have  poloidal  periods  1,  0  and  toroidal 
periods  0,  1  respectively  so  they  behave  like  the  angle  variables  in 
toroidal  geometry. 

Let  a   =   const,  be  a  drift  surface.   Then 


s  +  A„ 


where   A..   is   of   the  order  of  the  gyroradius  and  a   satisfies  the  first 
order  partial  differential  equation 


V  ,  •  V(7  =  0 

~d 


Here  the  drift  velocity  v,  is  given  by  (cf.  Section  IV. 1) 


Y^  -  ^  [B  ^  V  X  P||B] 


where  p..  -  Mv  /qB  is  the  parallel  gyroradius,  v  is  the  velocity 
component  of  the  particle  in  the  direction  of  the  magnetic  field,  M  is 
the  mass  and  q  is  the  charge.  The  form  for  the  drift  velocity  v^  has 
the   important  property  that  except  for  a  scalar  factor  it  is  divergence 
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free.  By  the  KAM  theory,  there  will  be  drift  surfaces  close  to  magnetic 
surfaces  if  the  gyroradius  is  sufficiently  small.  Substituting  the 
expression  for  v.,  neglecting  higher  order  terms  in  the  gyroradius,  and 
using  the  relations   B  •  Vs  -  J  •  Vs  -  0 ,  we  have 

B  •  VA,,  -  (B  X  7s)  •  Vp, 


Writing  VA,.  as 


dx,  ax.  dx. 

II    3s        dip  d<j> 


and  similarly  for  7p..  ,  we  see  that  X.    and  p.  satisfy  the  first  order 
equation 


(2.10) 


d  A II    d  p  n 


dxf)  dij> 


Here   we   have   used   the   two   representations   (2.8)  and  (2.9)  for  the 
magnetic  field  B.   From  the  relation 


„-Ilf''',W-.B.,*,l/^ 


between  the  energy  W  and  the  parallel  gyroradius  p.,  we  can  express  P||  in 
terms  of  the  electric  potential  *.  Here  fi  is  the  magnetic  moment  which 
is  an  adiabatic  invariant. 

1/2 
The  parallel  gyroradius  p   scales  like  M  ^  .   If  the  magnitude  of  * 

is  small  we  can  expand  p ,,  in  a  Taylor  series 


(2.11)  Pi,  =  A^  +  A,$  +  A„*^  + 
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Consider   the   term  i  -  i      (s)e^  '^in  the  Fourier  series  expansion 

mn  '^ 

for  the  electric  potential  *.   Substituting  this  into  (2.11)  and  making 
use  the  relation  (2.10),  we  obtain 


.     T   "•*  A,   imt^  -  i(n-4m)i/I_  , 
A  -  Ip   mn  1  e      ^  ^  ^'  P  + 

n-  tm 


This  exhibits  a  resonance  on  the  fliox  surface  s  -  Sg  where  the 
rotational  transform  is  the  rational  number  t  -  n/m.  To  find  the  island 
width  on  this  flux  surface  s  -  Sq ,  we  consider  the  modified 
representation  of  the  drift  surface  f(s  +  An )  -  const,  where  f  is  a 
function  to  be  chosen.  Since  s  is  of  order  one  and  X  is  a  small 
parameter  we  expand  f  in  a  Taylor  series 


(2.12)  f(s)  +  f'(s)A   +  ...  =  const. 


Since  A,,  has  a  singularity  at  s  =  Sq  ,  we  choose  f'(s)  -  s  -  s,,  in 
order  to  keep  the  second  term  of  (2.12)  finite.  Then  the  drift  surface 
has  the  representation 


2 
(2.13)  (s  -  Sq)  /2  +  (s  -  So)A..  =  const. 


Now  let 


emM 

cos  y 

n-  im      ^ 


where  «  is  a  scalar  factor  and  y  -  mV"  -  (n-  i.m)4>/l„   is  an  angle  variable. 
On  the  separatrix  s  =  Sq,  the  constant  on  the  right-hand  side  of  (2.13) 
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1/2 
is   eM  '    /t'(so)   by  I'Hospital's  rule.   Therefore   the  square  of  the 

island  width  is  given  by 


The   shear   t'Cs^)   should  not  be  zero,  which  recalls  the  conditions  in 

Moser's   theorem   (cf.   Section  I.l).   The  radial  step  size  scales  like 

1/4 
M  '    (cf.   [9]).    For  fixed  temperature  the  collision  time  step  scales 

1/2 
like  M  '    and  hence  the  transport,  which  scales  like  the  square  of  the 

step  size  over  the  time  step,  is  of  order  one.   This  means  that  at  equal 

ion  and  electron   temperatures,   the   loss   rates   are   of  comparable 

magnitude . 
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III.    STUDY  OF  MAGNETIC  SURFACES  ON  A  FINE  LENGTH  SCALE 


The  existence  of  nested  magnetic  surfaces  is  an  essential 
requirement  for  the  long  time  confinement  of  charged  particles  in 
toroidal  geometry.  The  existence  of  such  surfaces  has  only  been  proved 
rigorously  when  the  magnetic  field  has  special  2d  symmetry  [20].  In  the 
asjrmmetric  or  3d  case,  if  the  magnetic  field  has  shear  then  the  KAM 
theory  shows  that  most  of  the  magnetic  surfaces  persist  for  sufficiently 
small  and  smooth  perturbations  of  a  symmetric  case.  In  general,  the  KAM 
theory   only  applies  if  the  perturbation  is  extremely  small  and  there  is 

shear.   In  fact,  Arnold's  and  Moser's  proofs  require  the  perturbation  to 

-  333        -48 
be   less   than  10"    and  10    respectively  in  appropriate  units  [14]. 

Their  methods  are  vital  to  theoretical  dynamics,  but  are  numerically  of 

less  practical  value.   Numerical  evidence  suggests  that  the  perturbation 

can  be   sizable   [17],   but  an  upper  limit  is  not  yet  known.   From  the 

numerical   experiments,   surfaces   do   not   exist   in  some  models  if  the 

perturbation  is  large  enough  [17]. 

The  asymmetric  case  is  not  only  mathematically  interesting  but  is 
physically  more  realistic,  since  all  experimental  devices  are  to  some 
extent  asymmetric.  In  fact,  a  stellarator  is  by  its  very  nature 
asymmetric.  Since  a  small  break  in  symmetry  can  enhance  the  transport 
of  ions  and  electrons  significantly  it  becomes  important  to  study  the 
asymmetric  case  [8]. 

In  a  magnetic  field  the  motion  of  a  charged  particle  consists  of 
rapid  gyration  about  an  instantaneous  center,  called  the  guiding  center, 
which   closely   follows   the  magnetic  field  lines.   However,  the  guiding 
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center  orbits  drift  slowly  across  the  magnetic  field  because  of  its 
variation.  When  a  collision  occurs,  the  guiding  center  of  the  particle 
jumps  dlscontinuously  from  one  orbit  to  another  in  a  process  that 
generates  a  random  walk  across  the  field  lines.  The  rate  of  diffusion 
in  the  random  walk  is  proportional  to  the  square  of  the  step  size 
multiplied  by  the  collision  frequency.  For  good  surfaces,  the  step  size 
is  proportional  to  the  gyroradius  of  the  particle.  However,  in  the 
asymmetric  case  island  chains  or  stochastic  regions  may  form  in  the 
vicinity  of  rational  surfaces.  If  the  typical  width  w  of  the  island 
chain  or  the  stochastic  region  is  greater  than  the  gyroradius  of  the 
particle  and  there  are  very  many  islands,  then  w  becomes  the  step  size 
and  the  transport  is  enhanced.  Therefore  it  is  important  to  study  the 
regularity  of  the  magnetic  surfaces  and  drift  surfaces  on  the  length 
scale  of  the  gyroradius  of  the  particle. 

In  this  chapter  we  shall  investigate  these  issues  for  a  stellarator 
field,  in  particular  for  a  Heliac  field,  on  the  length  scale  of  the 
gyroradius  p  of  the  electron,  which  we  measure  in  units  of  the  plasma 
radius.  A  typical  laboratory  value  for  p  is  0.0002.  To  facilitate 
the  calculations  we  shall  use  an  especially  simple  analytical 
representation  of  the  magnetic  field  to  model  the  stellarator.  The 
solenoidal  character  of  the  field  will  be  preserved  everywhere.  The 
simple  representation  enables  us  to  trace  field  lines  around  the  torus 
many  times  with  limited  computing  resources.  In  the  following  study  of 
magnetic  surfaces,  we  usually  traverse  the  torus  thousands  of  times  in 
order  to  reveal  the  structure  of  the  magnetic  surfaces.  On  the  Cray 
computer,  Adam's  method  is  employed  to  integrate  the  field  lines  and  the 
integration   over   one  field  period  requires  only  30  /js .   A  typical  plot 
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of  magnetic  surfaces  or  ergodic  regions  such  as  appears  in  Figs.  11-14 
takes  about  1.5  hours  of  CPU  time. 

1.   Field  line  equations  for  a  simple  stellarator 

The   magnetic   field  lines  are  described  by  the  conservative  system 
of  ordinary  differential  equations  (ODE) 


(3.1)  ^  -  B(x) 

dt 


The  magnetic  field  B  in  a  current  free  region  (7  x  B  -  0)  can  be 
expressed  as  the  gradient  of  a  scalar  potential,  B  -  7^  .  Since  B  is 
divergence  free,  (j)  is  harmonic.  If  B  is  periodic  in  z  with  period 
2n/a  in  cylindrical  coordinates,  the  system  will  be  topologically 
toroidal  with  2ir/a  as  the  circumference  along  the  axis  of  the  torus. 
Hence  1/a  may  be  considered  as  the  major  radius  of  the  torus.  In  this 
case  a  convenient  representation  of  the  scalar  potential  is 


(3.2)         <^  -  BqZ  +  A-^e   +  X  I  tik   I^(akr)sin(i9-akz) 

i  k  ak 


where  I.  is  the  modified  Bessel  function  of  order  2  .  The  model  we 
choose  is  somewhat  simpler  than  the  usual  MHD  model  of  the  torus,  but  it 
is  physically  and  practically  meaningful.  The  model  provides  us  with  an 
analytical  representation  of  the  magnetic  field  for  which  we  can 
integrate  the  field  lines  with  great  accuracy  in  order  to  study  the 
magnetic  surfaces  on  a  very  fine  length  scale.  In  the  following  we 
shall  use  one  or  two  terms  under  the  summation  sign.  This  allows  us  to 
examine  the  effect  of  3d  geometry  on  the  magnetic  surfaces  efficiently. 
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The  corresponding  magnetic  field  B  -  (B^.  B^ ,  B^)   in  cylindrical 
coordinates  Is  given  by 


d<i) 


\  -  _  -  I  I  A.,  I'(akr)sin(i(?-Qkz) 


o     \    dd)        A  ,    1  V  T-   iA 

^S   -  -  —  -  -il  +  -11  ik  I.(akr)cos(i«-akz) 

r  36  r    r  i  k   ok 


^z  -  —  =  Bq  "  ^  ^  A   I  (Qkr)cos(id-akz) 


The  potential  d  alone  yields  the  field  B^  -  1/r  of  a  filamentary 
current  located  at  r  -  0.  The  introduction  of  the  potential  6  enables 
us  to  obtain  a  Heliac  stellarator  field  with  a  helical  excursion  of  the 
magnetic  axis. 

To  gain  insight,  let  us  first  consider  the  case  of  helical  symmetry 
where  B  depends  only  on  r  and  u  -  d-az  .  Translationally  and 
rotationally  invariant  examples  can  be  considered  as  limiting  cases. 
Here  magnetic  surfaces  can  be  found  by  introducing  a  flux  function  \j)  = 
0(r,u)  satisfying 

B  •  VV>  -  0 
Together  with  V  •  B  =  0,  this  gives 


^\f.^^\-^^\^f.-' 


-^ —  +  T'   (B.-  arB  )  -  0 

ar   au  ^  ^    z 
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Accordingly,  \&  satisfies  the  equations 

r  du      r 


ff-».-«^\ 


which  can  be  solved  to  yield 


0  A       ' 

V)  =  -A  ,log  r  +  B„Qr  -    r  Y^  _i   I  .(Qir)cos(iu) 

2       £      £      ^ 


Since  V  is  a  first  integral  of  the  ODE  system  (3.1),  the  problem  of 
finding  magnetic  surfaces  is  reduced  to  the  easier  question  of  finding 
level  curves  of  V  iri  two  -  dimensions . 

To  find  the  surfaces  of  V  =  constant  we  solve  the  ODE 


dr    . 

dt  =  "^u 


du      , 
d^  -  ■  ^r 


This  is  a  Hamiltonian  system  with  one  degree  of  freedom  that  is  always 
integrable.  Hence  the  existence  of  magnetic  surfaces  is  established  in 
this  case.  Fig.  1  shows  the  level  curves  for  i/i  =  constant  at  the  cross 
section  v  =-  0 .  Here  v  =  za/2n  is  the  normalized  toroidal  variable 
with  period  1.  The  magnetic  axis  and  the  separatrix  can  be  found  by 
solving 


i/.  =  0    ,       V  =  0 
^u  r 
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Fig.  2  shows  the  same  magnetic  surfaces  produced  by  a  line  tracing 
code  at  four  different  cross  sections  (cf .  Section  4) .  We  see  that  each 
cross  section  is  the  rotation  of  the  previous  cross  section  through  90 
in  the  counterclockwise  direction.  The  cross  section  at  v  -  0  of  Fig.  2 
is  identical  with  Fig.  1.  The  two-dimensional  case  can  also  serve  as  a 
check  on  the  line  tracing  code. 

2 .   Formation  of  stochastic  regions  in  three -dimensions 

If  we  add  a  harmonic  with  different  helicity,  for  example  a  term 
involving  A-,,  to  the  two-dimensional  potential  <f) ,  then  the  resulting 
system  is  three-dimensional  and  B  has  no  symmetry.  In  this  case  there 
is  no  flux  function  and  in  fact  magnetic  surfaces  do  not  in  general 
exist  everywhere.  Some  of  the  field  lines  may  sweep  out  complicated 
geometry,  such  as  island  chains  and  stochastic  regions. 

Figs.  3-7  show  a  sequence  in  the  development  of  a  stochastic  region 
near  the  separatrix  as  we  increase  the  magnitude  of  A„,  =  DEL21.  The 
stochastic  region  in  each  graph  is  formed  by  one  field  line.  For  A„, 
=  -.00001  we  cannot  see  the  stochastic  region  on  the  scale  of  the  plasma 
radius.  We  have  enlarged  a  portion  of  the  graph  near  the  separatrix  in 
Fig.  3b.  We  see  that  the  stochastic  region  prevails  there.  To  prove 
that  the  stochastic  region  is  not  due  to  numerical  errors  we  have  tried 
several  accuracies  for  the  ODE  solver  in  this  case,  and  the  results  are 
Che  same.  Therefore  we  conclude  that  the  surfaces  near  the  separatrix 
will  be  broken  into  stochastic  regions  no  matter  how  small  the 
perturbation  is.  The  stochastic  region  grows  out  from  the  separatrix 
exponentially  as  A„,  increases. 
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Fig.  8  shows  the  cross  sections  at  v  -  0  of  the  magnetic  surfaces 
for  A2,  -  -.11,  with  the  rotational  transform  t  increasing  from  0.25 
at  the  magnetic  axis  to  0.33  before  the  surfaces  break  up  at  the 
stochastic  region.  It  also  displays  an  island  chain  near  t  =  2/7. 
Outside  the  stochastic  region  it  appears  that  we  have  many  regular 
surfaces  again.  The  only  low  order  resonant  surface  corresponds  to  the 
rotational  transform  t  =  2/7.   This  is  a  case  we  shall  study  in  detail. 

We  find  that  a  large  stochastic  region  has  developed  near  the 
separatrix,  and  there  are  island  chains  where  small  stochastic  regions 
are  formed  in  the  region  nearer  the  magnetic  axis.  As  the  magnitude  of 
the  asymmetry  increases,  the  stochastic  region  near  the  separatrix  grows 
much  faster  than  the  width  of  the  islands  inside  the  region  around  the 
magnetic  axis.  We  conclude  that  we  have  good  flux  surfaces  close  to  the 
magnetic  axis,  but  most  of  the  surfaces  are  destroyed  near  the 
separatrix. 

In  the  literature  about  the  destruction  of  magnetic  surfaces  in  a 
stellarator  field,  the  field  lines  near  the  separatrix  in  the  perturbed 
system  usually  escape  to  infinity  [11].  However,  in  our  case  the  field 
lines  near  the  separatrix  stay  finite.  In  fact,  a  bounded  stochastic 
region  is  formed  by  following  the  orbit  of  a  single  field  line  near  the 
separatrix.  Outside  the  stochastic  region,  we  have  good  surfaces 
surrounding  it  again.  Therefore  the  stochastic  region  formed  by 
following  the  field  line  must  be  invariant.  We  also  see  that  there  are 
gaps  in  the  stochastic  region.  Thus  if  we  follow  the  field  line 
starting  from  another  point,  a  stochastic  region  disjoint  from  the 
previous  one  will  form  (cf.  Figs.  9a,  9b).  The  existence  of  such 
generalized  invariant  sets  has  only  been  proved  rigorously  by  Aubry  and 
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Mather  in  simple  cases  [22].  The  boundedness  of  the  stochastic  region 
we  have  exhibited  is  associated  with  the  simple  form  of  the  magnetic 
field  we  have  chosen  in  order  to  facilitate  numerical  integration  of  the 
ODE  system  (3.1). 

3 .   Width  of  islands 

As  we  observed  at  the  beginning  of  the  chapter,  the  formation  of 
island  chains  and  stochastic  regions  may  enhance  transport.  Therefore 
it  is  important  to  know  how  large  the  island  width  is. 

Islands  tend  to  form  in  the  vicinity  of  the  surfaces  with  rational 
rotational  transform  i  =  n/m.  To  find  the  width  of  an  island  we  have  to 
locate  its  center  first.  The  center  of  the  island  is  an  elliptic  fixed 
point  of  t",  the  m-fold  iterate  of  the  Poincare  mapping  T  which  is 
defined  by  following  the  field  lines  once  around  the  torus.  The 
elliptic  fixed  point  satisfies  the  equation  T  (x)  =•  x,  which  can  be 
solved  by  any  standard  method  in  numerical  analysis. 

After  locating  the  center  of  the  island,  we  pick  a  neighboring 
point  and  trace  from  it  a  magnetic  field  line  around  the  torus.  If  the 
magnetic  field  line  remains  close  to  the  fixed  point  after  every  m  turns 
around  the  torus,  then  this  field  line  will  form  an  island  around  the 
fixed  point.  We  can  repeat  the  procedure  of  tracing  the  magnetic  field 
line  for  points  further  out  until  a  stage  is  reached  where  the  magnetic 
field  line  fails  to  enclosed  the  fixed  point.  From  this  we  can  estimate 
the  width  of  the  island. 
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The  following  table  shows  the  island  width  w  in  units  of  p  for 
various  values  of  m.  Included  are  a  sample  of  the  low  order  rational 
surfaces.  We  find  that  the  island  width  decays  exponentially  with  m  -►  «. 


m 


w 


7 

10 

13 

17 

24 

31 

32.34 

15.20 

8.05 

1.42 

0.16 

0.03 

3.48 

2.72 

2.09 

0.35 

-1.83 

-3.51 

log  w 


Table  1.   Island  width  at  the  rational  surfaces. 

Fig.  10  establishes  the  exponential  decay  of  the  island  width. 
From  the  table  we  see  that  the  island  width  is  less  than  one  gyroradius 
of  the  electron  for  m  >  17 ,  so  that  higher  order  islands  are 
insignificantly  small  and  presumably  do  not  affect  transport.  The 
rational  surfaces  with  rotational  transform  t  =  n/m  whose  denominator  m 
is  less  than  17  are  not  numerous  in  the  range  0.25  <  t  <  0.33. 
Therefore  the  overall  effect  of  island  formation  on  transport  should  be 
very  small.  In  particular,  island  formation  due  to  the  3d  asymmetry  of 
the  magnetic  field  cannot  explain  the  anomalous  transport  of  electrons. 
We  shall  seek  for  other  explanations  in  the  next  chapter. 

The   island  width   at  a  rational  surface  with  rotational  transform 

n/m  has   been   shown  in  the  literature  to  be  proportional  to  the  square 

root   of   the  coefficient  a   of  a  Fourier  series  expansion  of  the  field 

mn 

with  respect  to  toroidal  and  poloidal  angles  [24].  Since  the  field  is 
analytic,  the  coefficient  in  the  Fourier  series  expansion  decays 
exponentially  with  m,  which  implies  the  exponential  decay  of  the  island 
width.    Our   calculation  confirms  this  theoretical  result  and  exhibits 
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quantitatively   the   rate   of  the  decay  in  a  typical  stellarator  example 
that  is  not  without  practical  interest. 

4.   Flux  surfaces  on  the  length  scale  of  the  electron  gyroradius 

To  substantiate  the  result  about  exponential  decay  of  the  island 
width  that  we  have  presented  in  Section  3,  we  use  the  line  tracing 
method  again  to  find  out  how  regular  the  magnetic  surfaces  are  on  the 
length  scale  of  the  gyroradius  of  the  electron.  We  choose  four  points 
that  are  one  electron  gyroradius  apart  in  the  plasma  region  and  follow 
the  field  lines  emanating  from  them  for  more  than  32000  cycles  around 
the  torus.  We  then  enlarge  the  intersection  of  their  orbits  in  the 
plane  v  -  0  to  bring  them  up  to  the  scale  of  the  gyroradius  of  the 
electron;  see  Figs.  11-14.  We  observe  that  if  they  are  not  close  to  a 
low  order  rational  surface,  the  surfaces  swept  out  behave  regularly. 

In  this  situation,  it  is  not  possible  to  predict  poor  confinement 
of  the  electrons,  since  the  stochastic  regions  are  very  small  and  are 
surrounded  by  good  magnetic  surfaces.  Recently  Bauer,  Betancourt  and 
Garabedian  have  developed  a  Monte  Carlo  code  to  calculate  neoclassical 
plasma  transport  [4,5].  The  numerical  results  suggest  that  anomalous 
transport  of  electrons  may  be  due  to  small  resonant  terras  in  the 
electric  potential.  The  Monte  Carlo  code  is  based  on  calculations  from 
a  three-dimensional  magnetohydrodynamic  equilibrium  code.  It  is  our 
goal  to  use  an  independent,  simple  analytical  formulation  to  study  this 
issue.   In  the  next  chapter  we  shall  discuss  this  matter  in  detail. 
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Figure  1.   Level  curves  of  the  flux  function  ^  -   constant  with 
ASPECT  -  2,  DEL-1  -  .1  and  DELll  -  .185,  where  ASPECT 
is  the  ratio  of  the  major  radius  to  the  plasma  radius. 
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Figure  2.   Magnetic  surfaces  produced  by  line  tracing 
at  V  -  .00,  .25,  .50,  .75  with  ASPECT  -  2.. 
DEL-1  -  .1  and  DELll  -  .185. 
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Figure  3a.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1,  DELll  -  .185  and  DEL21  -  -.00001. 
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Figure  3b.   An  enlargement  by  a  factor  of  200  of  Fig.  3a 
near  the  separatrix. 
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Figure  4.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1,  DELll  -  .185  and  DEL21  -  -.0001. 
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Figure  5.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1.  DELll  -  .185  and  DEL21  -  -.001. 
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Figure  6.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1,  DELll  -  .185  and  DEL21  -  -.01. 
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Figure  7.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1,  DELll  -  .185  and  DEL21  -  -.1. 
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Figure    8.      Magnetic    surfaces    of   a  Heliac    at  v  -   0  with   ASPECT   =   2, 

DEL-1   =    .1,    DELll   -    .185   and  DEL21   -    -.11.    The    stochastic 
region    is   bounded    inside   and  outside   by   good   surfaces. 
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Figure  9a.   Trajectory  of  one  field  line  with  ASPECT  -  2, 
DEL-1  -  .1,  DELll  -  .185,  DEL21  -  -.11  and 
initial  point  at  r  -  .85,  9-0. 
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Figure  9b.   Trajectory  of  one  field  line  with  ASPECT  -  2 
DEL-1  -  .1,  DELll  -  .185.  DEL21  -  -.11  and 
initial  point  at  r  -  .36,  6   -  it . 
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Figure  10.   Exponential  decay  of  the  island  width 
w  as  a  function  of  the  denominator  m 
of  the  rotational  transform  t    -   n/m. 
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Figure  11.   Magnetic  surfaces  of  k    field  lines  with  the 
distance  between  initial  points  equal  to  one 
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Figure  12.   Magnetic  surfaces  of  4  field  lines  with 
negative  slope. 
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Figure  13.   Magnetic  surfaces  of  4  field  lines  with 
positive  slope. 
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Figure  1^.   Magnetic  surfaces  of  4  field  lines  with 

initial  points  one  electron  gyroradius  apart. 
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IV.   DRIFT  SURFACE  ISLANDS 


The  results  in  the  previous  chapter  suggest  that  the  formation  of 
island  chains  due  to  the  asynunetry  of  the  magnetic  field  cannot  explain 
the  anomalous  transport  of  electrons  observed  in  experiments.  In  this 
chapter,  we  continue  to  investigate  this  issue  by  computing  the  guiding 
center  orbits  of  the  electrons  and  estimating  island  widths  of  the  drift 
surfaces  that  are  swept  out.  The  effect  of  the  electric  field  is 
included  in  the  guiding  center  equations.  We  suspect  that  the  anomalous 
transport  of  electrons  is  due  to  resonance  of  the  electric  field.  We 
shall  study  how  the  electric  potential  affects  the  island  width 
numerically  to  substantiate  this  conjecture. 

At  this  stage  it  is  appropriate  to  give  a  brief  description  of  the 
guiding  center  theory  of  a  charged  particle  in  an  electromagnetic  field 
in  order  to  obtain  a  better  understanding  of  the  results  presented  here. 
In  the  following  we  shall  assume  that  magnetohydrodynamic  equilibrium 
has  been  reached  so  that  the  fields  are  static. 

1 .   Guiding  center  equations 

The  equation  of  motion  of  a  charged  particle  in  a  magnetic  field  B 
and  an  electric  field   E   is  given  by 


M^-q  (vxB+E) 
dt 
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Here  v  is  the  velocity  of  the  particle,  M  is  the  mass  and  q  is  the 
charge.  When  the  magnetic  field  and  the  electric  field  are  uniform  in 
space,  the  equation  of  motion  can  be  integrated  explicitly.  The  motion 
of  the  charged  particle  is  found  to  be  a  superposition  of  a  rapid 
circular  gyration  around  an  instantaneous  center,  called  the  guiding 
center,  with  cyclotron  frequency  u>  -  qB/M,  and  a  drift  motion  across  the 
magnetic  field  lines.  Here  B  is  the  magnitude  of  B.  The  radius  of 
gyration  is  called  the  gyroradius . 

When  the  magnetic  field  and  the  electric  field  are  nonuniform 
because,  for  example,  the  magnetic  field  has  a  spatial  gradient  so  the 
magnetic  field  lines  are  curved,  the  integration  of  the  equation  of 
motion  cannot  in  general  be  carried  out  in  closed  form.  However  for  most 
cases  of  interest  in  fusion  research,  the  scale  of  nonuniformity  of  the 
fields  is  very  large  compared  to  the  gyroradius,  and  the  time  scale  is 
very  long  compared  to  the  gyroperiod  l-n/w.  Under  these  conditions  we 
can  expand  the  B  field  and  the  E  field  about  the  guiding  center. 
Writing  the  motion  of  the  charged  particle  as  the  sum  of  the  gyration 
and  the  motion  of  the  guiding  center  and  taking  a  time  average  over  one 
gyroperiod,  we  obtain  the  guiding  center  equations  in  the  form  [19] 


(4.1)  ^-  -  Zl    [(1-P  k*Vxb)B  +  V  X  p  B]  +0(e^) 

dt    B      "  " 


where  x  is  the  position  vector  of  the  guiding  center,  v..  is  the 
component  of  the  velocity  in  the  direction  of  B,  P ,,  ="  Mv  /qB  is  the 
parallel  gyroradius,  and  b  is  the  unit  vector  in  the  direction  of  B. 
Here  e    is  the  ratio  of  the  gyroradius  to  the  average  plasma  radius. 
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Note   that   Piib'Vxb   and  Vxp.,B   are   of  the  order  t.      Dividing  the 
guiding  center  equations  by  the  factor  (l-P||b«Vxb)  ,  we  have 

1    dx    V 


^  -  IJL  [  B  +  ^  ^  ^11^   ] 


1-p  b'Vxb  dt    B        l-P||b«Vx^ 


Expanding   the  right-hand  side  of  the  resulting  equations  in  powers  of  « 

2 
and  neglecting  terms  of  order  €  ,  we  obtain 


(4.2)  ^-Y-;;;i[B  +  Vx  p,.B] 

dt    ^   B  " 


Here   the   time  variable  has  been  rescaled,  but  the  trajectory  of  the 

guiding  center   does  not  change  because  the  right-hand  side  of  (3.2)  is 

2 
independent  of  time.    The  new  equation  is  also  accurate  to  order  t    . 

Therefore  we  do  not  lose  anything  by  transforming  (4.1)  to  (4.2). 


The  form  for  the  drift  velocity  v  in  (4.2)  has  an  important 
property  related  to  the  KAM  theory.  Except  for  a  scalar  factor,  which 
is  of  no  consequence,  it  is  divergence  free.  Since  the  ordinary 
differential  equations  (4.2)  are  autonomous,  the  guiding  center  orbit 
has  the  same  trajectory  as  the  divergence  free  vector  field  B  +  7  x  p,,^- 
This  means  that  Poincare  map  defined  by  (4.2)  is  measure-preserving  and 
the  KAM  theory  applies.  The  vector  field  B  +  V  x  p  B  is  a  perturbation 
of  the  magnetic  field  B  and  the  perturbation  V  x  p..B  is  of  the  order  of 
the  gyroradius.  The  KAM  theory  suggests  that  if  the  gyroradius  is 
sufficiently  small  and  magnetic  surfaces  exist,  then  there  are  nearby 
drift  surfaces. 

Throughout  the  motion  the  energy 
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1   2 
W  =  T  Mv,.  +  /iB  +  q* 


1    2 
is   conserved  and   the   magnetic  moment  /i  -  -  Mv  /B   is  an  adiabatic 

invariant.  Here  v  is  the  component  of  y  perpendicular  to  B  and  * 
is  an  ambipolar  electric  potential,  E  -  V*.  The  constancy  of  the  energy 
W  is  a  consequence  of  the  equation  of  motion  and  the  fact  that  the 
electric  field  is  static.  The  invariance  of  the  magnetic  moment  /i 
follows  from  the  condition  that  the  change  of  the  magnetic  field  in  a 
gyroradius  is  small  compared  with  the  magnetic  field  itself.  The  guiding 
center  equations  are  obtained  by  the  method  of  averaging  over  the 
gyroperiod.  The  method  of  averaging  leads  directly  to  the  calculation 
of  all  adiabatic  invariants  which  are  approximate  integrals  of  the 
motion  obtained  by  averaging  over  the  fast  gyration,  and  one  can  expand 
in  a  similar  way  to  analyze  bounce  motion  and  then  drift  motion  (17]. 
The  adiabatic  invariants  corresponding  to  the  bounce  motion  and  the 
drift  motion  are  the  longitudinal  invariant  and  the  flux  invariant 
respectively.  However,  we  shall  deal  directly  with  numerical 
integration  of  the  guiding  center  equations. 

The   constancy  of  the  magnetic  moment  and  the  energy  will  determine 
v  as  a  function  of  position  given  by 


2  1/2 

V||  =  [^(W  -  mB  -  q*)]^ 


In  order  to  uniformize  the  branching  where  v  vanishes  at  turning 
points,  the  system  is  augmented  by  an  ordinary  differential  equation  for 
the  parallel  gyroradius, 
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(4.3)  ^  -  5  [  B-V  ipf  +  p,|V  X  B.V  Iph 

dt    M       2  "     "        2 


which  follows  from 


at 


1  2 
Here  Che  vector  V  —p        is  found  by  differentiating  W  to  obtain 


2         ^     q^B^         qB^ 


The  conservation  of  energy  can  be  used  to  test  the  accuracy  of  the 
numerical  solution.  The  magnetic  moment  /i  is  held  constant  throughout 
the  motion. 

It  is  convenient  to  write  equations  (4.2),  (4.3)  in  terms  of  non- 
dimensional  quantities  and  to  integrate  the  dimensionless  form  of  the 
differential  equations.  All  physical  quantities  are  scaled  according  to 
familiar  formulae  given  in  next  chapter.  After  substituting  all  the 
scaled  quantities  into  the  guiding  center  equations,  we  obtain  Che  orbit 
equations  in  a  dimensionless  form.  In  our  model  V  x  B  =  0,  so  the  above 
equations  reduce  to 


^-  P„B.  ((p.)2  -  (  4+^  )VB  -  !!  )  X  B 
dt     "  B    ^2       ^2 


^  =  B  .  ((p,)2  -  (  4+^  )VB  -  !!  ) 
dC  ^       B    ^2       g2 
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where  p,  is  the  ratio  of  the  gyroradius  to  the  average  plasma  radius. 
These  equations  are  simple  and  the  right  hand  side  of  the  equations  is 
easy  to  compute  for  the  elementary  form  of  B  that  was  described  in 
Chapter  III.  The  CPU  time  required  to  integrate  over  one  field  period 
is  50i/s  on  the  Cray  X-MP/48  computer. 

2 .   Drift  surfaces 

In  this  section  we  study  the  drift  surfaces  for  electrons  and  ions. 
The  effect  of  the  electric  potential  is  included  in  the  computations. 

We  have  followed  guiding  center  orbits  for  different  values  of  p  . 

Li 

For  ions  a  typical  value  of  p^  is  0.01.  Without  an  electric  potential 
the  drift  surfaces  behave  relatively  well.  With  a  small  electric 
potential,  which  is  only  a  few  percent  of  the  temperature,  we  find  that 
some  drift  surfaces,  especially  near  the  separatrix,  become  destroyed 
and  form  stochastic  regions.  This  means  that  particles  may  escape  more 
rapidly  in  the  presence  of  an  electric  potential. 

For  electrons  at  a  comparable  temperature,  p.  -  0.0002.  The 
electron  drift  surfaces  are  therefore  closer  to  the  magnetic  surfaces. 
When  there  is  no  electric  field,  they  are  quite  regular.  With  a  small 
electric  potential  most  of  the  drift  surfaces  remain,  but  the  width  of 
the  associated  islands  at  rational  surfaces  is  of  the  order  /p~  (cf. 
Section  II. 4).  The  existing  analytic  theory  does  not  fully  explain  the 
anomalous  electron  transport,  which  is  one  or  two  orders  higher  than 
what  is  usually  predicted  [16].  It  is  our  goal  to  give  a  plausible 
explanation  of  this  phenomenon  by  studying  the  enhancement  of  the  island 
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width  due   to   an  electric  potential.    We  will  show  how  the  electric 
potential  affects  the  island  width  in  the  next  section. 

3.   Island  width  and  the  electric  potential 

We  study  the  fine  structure  of  the  drift  surfaces  swept  out  by  the 
electron  guiding  center  orbits.  Using  the  same  algorithm  as  in  the  case 
of  the  magnetic  field,  we  can  compute  the  drift  surface  island  widths  at 
the  rational  surfaces.  In  this  study  island  widths  near  a  sample  of  the 
rational  surfaces  with  rotational  transforms  2/7,  3/10,  3/17,  5/23  and 
7/31,  which  are  all  between  0.25  and  0.33,  are  computed. 

With  no  electric  potential,  the  island  width  decays  exponentially 
with  m,  the  denominator  of  the  rotational  transform  t  -  n/m  of  the 
rational  surface.  The  decay  rate  is  more  or  less  the  same  as  in  the 
case  of  the  magnetic  field.  In  fact  the  island  width  is  less  than  one 
gyroradius  of  the  electron  for  m  greater  17.  Therefore  higher  order 
islands  do  not  affect  the  transport  significantly. 

It  is  customary  to  find  the  electric  potential  *  from  Ohm's  law.  A 
conventional  argument  suggests  that  the  electric  potential  depends  on 
the  flux  variable  s  only.   Indeed,  consider  the  Ohm's  law 


(4.4)  E  +  yxB-rjJ 


where   n    is   the   resistivity  and   J    is   the  current  density.   The 

'o  ■' 

resistivity  r\      is  in  general  very  small  in  fusion  applications.   Taking 
a  scalar  product  with  B,  (4.4)  becomes 
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E  •  B  -  T?  J  •  B 
-   -    'o    - 


which  means  that  the  parallel  electric  field  is  the  product  of  the 
resistivity  and  the  parallel  current  density.  If  the  resistivity  rj  is 
very  small  or  vanishes,  then  E  —  74  is  perpendicular  to  the  magnetic 
flux  surfaces.  This  means  that  *  can  only  depend  on  the  flux  variable  s 
when  I  is  irrational.  However,  if  we  consider  the  generalized  Ohm's  law 
proposed  by  Boozer  [7] 


fi        J.B 

(4.5)  E  +  V  X  B  -  r,oJ  -  -  V.(r,  7^) 

B      ^   B 


where  r;.  is  a  scalar,  the  second  term  on  the  right  hand  side  of  (4.5), 
which  models  fluctuations  of  the  magnetic  field,  is  not  in  general  so 
small  (cf.  [13]).  In  this  case  the  above  argument  fails  to  assert  that 
the  electric  potential  *  is  a  function  of  the  flux  variable  only. 

Taking   a   scalar  product   of  equation  (4.5)  with  B  and  using  the 
divergence  theorem 


J.B 

J  V.(»,^V— )  dx  -  0 


we  have 

(4.6)  J"  ^•(<1'B)  dx  -  J  E'B  dx  -  /  "l^-^'^   dx 

which  means  that  the  loop  voltage  is  the  product  of  the  resistivity  and 
the  net  current  after  (4.6)  is  divided  by  the  flux.  This  enables  one  to 
calculate   the   resistivity  rj        by  measuring   the  loop  voltage  and  the 
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current  experimentally.  Note  that  (4.4)  also  implies  (4.6).  This  means 
that  with  the  generalized  Ohm's  law,  one  has  more  freedom  to  specify  the 
form  of  the  electric  potential  without  altering  the  relationship  between 
say  loop  voltage  and  net  current. 

In  our  formulation,  cylindrical  coordinates  are  employed.  It  is 
difficult  to  express  the  physical  quantities  numerically  in  terms  of 
flux  coordinates.  To  simplify  the  computation  we  introduce  an  electric 
potential  of  the  form 


*  -  *.,T  cos(i9-kz) 


in  cylindrical  coordinates,  where  T  is  the  temperature,  i  and  k  are 
integers  and  *.,  is  some  amplitude.  Since  the  electric  potential  does 
not  depend  only  on  the  flux  variable,  the  above  form  for  the  electric 
potential  serves  our  purpose  well,  as  we  shall  see  the  computation 
below. 

We  have  computed  the  island  widths  of  drift  surfaces  near  the 
rational  surfaces  for  several  different  values  of  *.,  with  i  -  17  and 
k  =-  10.   The  results  are  given  in  Table  2. 
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11*11  "  a 

7       10      17     23     31 

0.0000  21.88  10.16  1.41  1.06  0.03  0.250 

0.0025  25.00  13.10  5.78  2.63  0.25  0.180 

0.0100  20.94  14.06  11.25  4.45  0.81  0.128 

0.0200  20.20  16.90  15.30  6.00  1.60  0.103 

0.0400  20.31  26.50  23.80  6.50  3.50  0.084 

0.0600  20.00  25.50  23.80  5.50  5.50  0.069 

Table  2.   Drift  surface  island  v;idth  w  and  the  decay  rate  a   at 
different  values  of  | |*| | . 

The  island  width  is  measured  in  units  of  the  gyroradius  of  the  electron 
and  the  magnitude  of  the  electric  potential  |  |*|  |  -  l*pi.l  is  measured  In 
units  of  the  temperature  T,. 

For  each  value  ||*||,  the  island  width  obeys  the  exponential  law 
w  =■  Ae"'*™  for  some  decay  rate  q  (cf.  Figs.  15-18).  We  find  that  1/a 
varies   linearly  with  respect  to  the  square  root  of  the  magnitude  of  the 


electric  potential  J\ |*| | ;  see  Fig.  19.  This  means  that  the  island 
width  decays  exponentially  at  a  slower  rate  when  the  magnitude  of  the 
electric  potential  gets  larger. 

Since,  for  a  large  collection  of  islands,  the  total  width  is 
proportional  to  1/a,  the  step  size  is  approximately  proportional  to  the 
square  root  of  the  magnitude  of  the  electric  potential,  which  represents 
the  perturbation  of  the  system.  For  fixed  temperature  both  the  magnitude 
of  the  perturbation  5A  of  the  drift  surfaces  associated  with  $  and  the 
collision   time   scale   like   M  ^    ,    the  square  root  of  the  mass  M  of  the 
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particle.  Therefore  the  random  walk  step  size  scales  like  M  .  Hence 
the  transport,  which  scales  like  the  square  of  the  step  size  over  the 
time  scale,  is  of  order  one,  independent  of  the  mass  of  the  particle. 
This  means  that  with  equal  temperatures  the  transport  for  ions  and 
electrons  is  of  the  same  order  of  magnitude  and  the  condition  for  charge 
neutrality  may  be  more  readily  satisfied. 

Note  that  the  island  width  for  the  m  -  7  drift  surface  is  sizable 
when  the  electric  potential  $  is  zero.  This  means  that  the  formation  of 
the  islands  for  this  drift  surface  is  due  to  the  asymmetry  of  the 
magnetic  field.  The  electric  field  is  just  a  small  perturbation  to  the 
system,  so  that  the  island  width  does  not  change  much  at  this  drift 
surface.  However,  the  island  widths  change  appreciably  at  the  m  =-  17 
drift  surface  and  other  drift  surfaces  as  the  magnitude  of  $  increases. 
This  means  that  the  electric  potential  has  enhanced  the  higher  order 
island  widths .  Resonance  may  occur  there  and  this  may  cause  anomalous 
transport. 

In  our  numerical  experiments,  the  magnitude  of  the  electric 
potential  is  only  a  few  percent  of  the  temperature,  which  could  not  be 
detected  in  an  actual  experiment.  The  above  results  suggest  that  a 
small  electric  potential  may  cause  the  breakdown  of  drift  surfaces  and 
enhance  transport.  In  the  next  chapter  we  shall  use  the  Monte  Carlo 
method  to  calculate  the  confinement  time  of  the  electrons  and  compare 
results  with  and  without  the  electric  potential. 
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Figure  15.   Exponential  decay  of  the  island  width 
when  the  electric  potential  *  -  0. 
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Figure  16.   Exponential  decay  of  the  island  width  when 

the  electric  potential  *  -  0 .01Tcos(17^-10z) 
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Figure  17.   Exponential  decay  of  the  island  width  when 
the  electric  potential  *  -  0.02Tcos(17tf -lOz) 
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Figure  18.   Exponential  decay  of  the  island  width  when 

the  electric  potential  *  -  0.04Tcos(17tf-10z) 
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Figure  19.   Proof  that  the  reciprocal  of  the  decay  rate  a 
varies  linearly  with  the  square  root  of  the 
magnitude  of  the  electric  potential  *. 
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V.   MONTE  CARLO  CALCULATION  OF  CONFINEMENT  TIME 

In  this  chapter,  a  Monte  Carlo  code  to  calculate  neoclassical 
transport  for  ions  and  electrons  is  developed.  The  goal  of  the  code  is 
to  substantiate  the  results  we  have  obtained  in  Chapter  IV.  The  main 
conclusion  is  that  introduction  of  a  small  electric  potential  may 
increase  the  island  width  and  enhance  the  transport  of  electrons 
significantly. 

In  the  Monte  Carlo  calculation  of  neoclassical  transport,  a  back- 
ground equilibrium  magnetic  field  is  assumed  to  be  given.  Our  code  is 
developed  in  the  same  spirit  as  the  code  BETAMOC  by  Bauer,  Betancourt 
and  Garabedian  [4] .  In  BETAMOC  the  magnetic  field  is  obtained  from  the 
three-dimensional  magnetohydrodynamic  equilibrium  code  BETA  [3].  In 
our  code,  the  magnetic  field  is  calculated  more  directly  by  solving 
V  X  B  -  0  and  V  •  B  =  0.  A  simple  analytical  representation  of  B  which 
is  easy  to  compute  in  cylindrical  coordinates  has  been  described  in 
Chapter  III. 

The  Monte  Carlo  algorithm  was  first  described  by  Boozer  and  Kuo- 
Petravic  [8].  In  their  formulation  it  is  assumed  that  a  nested  family 
of  magnetic  flux  surfaces  exists.  As  in  the  BETA  code,  the  magnetic 
field  B  is  represented  as  B  =  Vs  x  Vrp,  where  s  and  0  are  flux  functions. 
The  function  s  is  assumed  to  be  single-valued  and  to  define  a  nested 
family  of  toroidal  surfaces  s  =  constant  for  values  of  s  ranging  from 
s  =  0  at  the  magnetic  axis  to  s  -  1  at  the  outermost  magnetic  surface. 
In  our  model  all  physical  quantities  are  computed  in  cylindrical 
coordinates,   but  we   have   developed  a  method  to  find  the  radial  flux 
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variable  s,  which  is  important  in  Monte  Carlo  calculations.  The  method 
uses  knowledge  of  the  position  of  the  magnetic  axis  and  magnetic 
surfaces  and  will  be  described  in  Section  4. 

In  the  following  section  we  shall  describe  the  Monte  Carlo 
algorithm  to  solve  the  drift  kinetic  equation  from  which  confinement 
time  and  transport  coefficients  can  be  estimated.  The  method  is 
basically  that  of  Boozer  and  Kuo-Petravic  [8]  and  Bauer,  Betancourt  and 
Garabedian  [5],  The  main  difference  lies  in  our  representation  of  the  B 
field. 

1 .   The  Monte  Carlo  algorithm 

In  the  drift  kinetic  theory  the  motion  of  the  particles  is 
described  by  the  Fokker- Planck  equation 


(5.1)  ^t  "■  ^d  •  ^^  -  2  •'d^f^l 


where  f  -  f(x,v,t)  is  the  distribution  function  for  particles,  u,  is  the 
collision  frequency,  y,  is  the  guiding  center  drift  velocity,  and  y  is 
the  particle  velocity.  Here  7  stands  for  gradient  in  physical  space  and 
L  is  a  second-order  partial  differential  operator  in  velocity  space  that 
approximates  the  collision  operator.  The  full  collision  operator 
changes  both  the  energy  and  the  magnetic  moment  of  the  particles.  In 
neoclassical  transport  collisions  are  dominated  by  pitch  angle 
scattering  associated  with  a  Lorentz  collision  operator  represented  by 
the  Laplacian  A  on  the  unit  sphere   in  velocity  space.    In  our 
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computations   energy  scattering   is   also   included.    We  shall  discuss 
details  about  the  collision  operator  in  Section  2. 

The  basic  idea  in  solving  the  Fokker-Planck  equation  is  as  follows. 
Initially,  a  large  collection  of  particles  are  placed  in  physical  space 
with  a  prescribed  distribution  of  energy  and  magnetic  moment.  The 
particles  are  allowed  to  move  along  the  guiding  center  orbits  with 
energy  and  magnetic  moment  held  constant.  After  a  time  step  related  to 
the  collision  frequency,  the  particles  are  subjected  to  collisions  which 
change  the  energy  and  magnetic  moment  according  to  rules  that 
approximate  the  collision  operator.  Following  the  particles  for  many 
such  time  steps,  we  obtain  an  approximate  solution  to  (5.1).  By 
estimating  the  decay  of  suitable  functionals  of  f  we  can  compute  the 
particle  confinement  time  and  transport  coefficient. 

The  Fokker-Planck  equation  is  solved  in  two  stages  in  our  iteration 
scheme.  In  the  first  step,  the  solution  is  advanced  by  integrating 
along  the  characteristics  of  the  first  order  partial  differential 
equation 


(5.2)  f^  +  v,.Vf  =  0 

t   ~a 


which   represent   guiding   center   orbits.    In   the   second   step,   the 
particles  undergo  a  random  walk  that  models  the  diffusion  equation 


(5.3)  f^  -  i  u^L[f] 
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in  velocity  space.    Combining   the   two   steps   repeatedly  we  get  an 
approximate  solution  of  (5.1). 

The   solution   to   (5.2)  is  equivalent  to  following  particles  along 
the  characteristics,  which  are  the  drift  orbits 


dx 

dt  "  -d 


since  f  is  constant  along  characteristics.  The  solution  to  the  orbit 
equations  has  been  discussed  in  Chap-er  IV.  The  solution  of  (5.3)  is 
solved  by  a  random  walk  method  which  is  best  illustrated  by  a  simple 
example . 

Suppose  we  want  to  solve  the  one -dimensional  diffusion  equation 


(5.4)  u-Du     ,     -<»<x<<» 

C       XX 


u(x,0)  -  6(x) 

where   D   is  a  constant  and  5  is  the  Dtrac  delta  function.   The  solution 
is  the  Gaussian  density  function 


1  ^2 

(5.5)  u(x,t)  -  r—z   exp(-  t^) 

(ATTDt)^/^  ^ 


with  mean  zero  and  variance  2Dt.   Discretizing  the  diffusion  equation  by 
a  finite  difference  scheme,  we  have 
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n+1  n     n    «  n  .   n 

(5.6)  ^   -^     A+r  ^^  -^  ^-1 

(Ax) 


2 
If  (Ax)   -  2DAt  then  (5.6)  becomes 


n+1    In    ,  1   n 


which   describes  a  random  walk  on  the  x-axis.   Therefore  the  solution  of 

(5.4)  can  be  constructed  as  follows.  At  time  t  -  0 ,  we  place  N 
particles  with  equal  mass  density  1/N  at  x  -  0.  Each  particle  is  then 
allowed  to  perform  a  random  walk  on  the  x-axis  with  step  size  y2DAt . 
After  n  steps  the  position  of  each  particle  is  just  the  sum  of  its 
displacements.  For  large  N  the  mass  density  distribution  on  the  x-axis 
has  mean  zero  and  variance  n«2DAt  =-  2Dt  .  By  the  central  limiting 
theorem  the  mass  distribution  will  converge  to  the  Gaussian  distribution 

(5.5)  as  n  -►  «. 

For  general  initial  data  u(x,0)  =  f(x)  the  solution  can  be 
constructed  by  placing  N  particles  appropriately  along  the  x-axis.  The 
mass  distribution  along  the  x-axis  at  t  =  nAt  will  approximate  u(x,t) 
for  large  N  and  small  At. 

2 .   Random  walk  simulation  of  the  collision  operator 

In  diraensionless  form  the  general  drift  kinetic  equation  including 
pitch  angle  scattering  and  energy  scattering  is 


(5-7)  f,  -  v^-Vf  -  \u^M   .  -\   f^[v\(vf  .  I   f^) 

V 
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where  A  Is  the  Laplacian  on  the  sphere  in  velocity  space,  i/  and  i/  are 
collision  frequencies,  and  T  is  the  temperature.  The  variation  of  f  due 
to  pitch  angle  scattering  is  represented  by 


/c  ON  if   1  A<r  1   a_.,  2,af 


V  ri 

where   r;   -  _!!.   is   the  cosine  of  the  pitch  angle  and  f  is  assumed  to  be 

V 

independent   of   the   azimuthal   angle   in  velocity  space.   The  Lorentz 
collision  operator  is  approximated  by  the  random  walk 


(5.9)  r,^^^^  -  a-i.^St)v^  ±    [(l-fll)i^^St]^/^ 


Here  the  choices  of  the  plus  or  minus  sign  are  equally  distributed  by 
means  of  a  random  number  generator.  The  magnetic  moment  is  then  updated 
by  the  formula 


W      2      * 
M--(l  -  r,^)(l  -  J) 


To  derive  (5.9),  consider  a  random  walk  with  the  step  sizes 


Srj     -  est  ±  D/St 


where  5t  is  the  time  step  and  C  and  D  are  coefficients  to  be  determined. 
Write 


(5.10)  f(t+6t,rj)  -  ^  f(t,r;+5f7^)  +  ^  f(t,f7+6»7_) 
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which  means  that  the  distribution  function  f  at  r;  has  equal 
contributions  from  the  two  neighboring  points  at  the  previous  time  step. 
From  a  Taylor  expansion  (5.10)  becomes 


(6t)fj.  =  j(    Sr,^   +  5r,_)f^   +  J[(«'?+)   +  (SO  ]f^^  + 


Dividing   through  by  St,    comparing  with  the  coefficients  of  the  partial 

3/2 
differential   equation   (5.8),   and  neglecting  terms  of  order  (St)  or 

2       2 
higher,  we  find  that  C  -  -u.rj   and  D  =  (1-f;  )i/,,  which  implies  (5.9). 


For  energy  scattering  the  collision  operator  is  represented  by 


/=  -.ix  M   i_  a_,  2  ,  .  ^  T  if, 

(5-11)  J^  -  ~i^[^  ^e^^^  -^  S  av) 

V 


which  has  a  steady  state  Maxwellian  solution 


97rT   '  19 

f  -  (^)    exp(-  ^Mv^/T) 


1   2 
Let  E  -  xMv  be  a  new  independent  variable  so  the  equation  (5.11)  can  be 

written  as 


Then  the  rule  for  energy  scattering  becomes  [8] 


■\  F  ^''p  1/2 

(5.13)       E^^^^  =  E^    -  2.^6t[E^  '  (f  +  iT  ^^^^  ^  2(VeT^^) 

e 
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To   derive   (5.13)   we  need  a  more  general  type  of  random  walk.   As 
before  let  6E     -   C5t  +  UjSt,    but  write 


+ 


f(t+5t,E)  -  Pj^f(t,»;+5E^)  +  q^f(t,r}+SE_)    + 

P2f(t.»7+25E^)  +  q2f(t,r;+26E_)  +  rf(t,E) 

where   p..+q,+p^+q^+r-l .   Here  p,  and  q,  are  of  order  1,  p„  and  q^  are  of 

1/2 
order  (6t)    ,  and  r  is  of  order  St.      Using  a  Taylor  expansion  again,  we 

have 

(5t)f^  -  (p^6E^  +  qj^5E_  +  P225E^  +  q225E_)fg  + 


^[p^(6E^)^  +  q^(5E_)^  +  P2(26E^)^  +  q2(26E_)^]  f^g.  +  rf(t,E) 


Dividing  the  equation  by  5t  and  comparing  with  the  coefficient  of  f^p  in 

2  1       1/2 

(5.12),   we   get   D   -  4i/^ET.   If  we  let  P^  -  q^^  "  2  (l-"  (-^t)  '^  )  ,   q2  -  0 

1     1/2 
and   p^   -  TD(5t)    /T  then  by  comparing  with  the  coefficient  of  fp  we 

get  (5.13). 

3 .   Estimation  of  the  loss  rate 

In  this  section  we  present  a  method  to  estimate  the  particle 
confinement  time  given  by  Garabedian  [5].  The  problem  of  estimating  the 
confinement  time  or  loss  rate  becomes  difficult  because  of  the  wide 
range  of  collisionality  that  occurs  in  practice. 
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In  practice,  the  collision  frequency  is  not  high.  A  general 
concept  valid  for  the  full  range  of  collision  frequency  is  that  of 
exponential  decay.  The  idea  is  to  estimate  confinement  time  from  the 
exponential  decay  of  the  distribution  function  f.  Consider  the  mixed 
initial  and  boundary  value  problem  for  the  diffusion  equation 


Uj.  -  |^(Ds|j)         u(s,0)  -  F^(s).     u(l.t)  -  0 


with  regularity  imposed  as  a  boundary  condition  at  s  -  0.  For  this 
there  exists  a  complete  set  of  eigenfunctions  u  (s)  with  eigenvalues  A  . 
The  Green's  function  is 


-A  t 

G(s,(7,t)  -  Z  u  (a)u  (s)e  " 
n    n 


The  expected  value  for  any  function  F_  has  the  form 


1  N 

^  Z  F2(s  )  -  //  F^(a)F^(s)G(s,a.t)    dsda 

.  Ae   °   =  Ae-'^/^ 


where  A  is  a  scalar  factor.  This  will  provide  an  effective  estimate  of 
the  confinement  time  t  if  F,  and  Fj  are  good  approximations  to  the 
lowest  eigenfunction  u  of  the  problem.  In  our  case  we  use  ^i  ~  ^2  " 
1-s.   This  leads  to  the  expected  value 


H(t)  =  ^  S  (1-s  )  =  Ae"^/' 
1     ■' 
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where  s.  measures  the  radial  position  of  the  i th  particle  at  the  time  t. 
J 

We  compute  H(t)  at  each  time  t  and  fit  a  line  to  the  function  log  H(t) 
by  least  squares  to  obtain  the  slope  which  gives  an  estimate  of  the 
confinement  time  r. 

4.   A  method  of  evaluating  the  radial  position  of  a  particle 

In  the  Monte  Carlo  method  that  has  been  described  it  is  necessary 
to  find  the  radial  position  s  of  the  particles  in  order  to  calculate  the 
confinement  time.  In  our  model  a  simple  analytical  representation  of  the 
magnetic  field  in  cylindrical  coordinates  is  used  for  the  computations. 
Determination  of  the  radial  position  s  of  a  particle  is  not  that  direct 
as  in  the  BETAMOC  code,  where  s  is  one  of  the  variables.  We  use  a  line 
tracing  method  to  find  the  magnetic  axis  and  the  outermost  magnetic 
surface,  from  which  an  approximation  for  the  radial  coordinate  s  can  be 
obtained.   We  give  here  a  brief  description  of  the  algorithm  to  find  s. 

We  first  find  the  magnetic  axis,  which  is  a  field  line  closing  on 
itself  after  one  field  period,  and  we  locate  the  outermost  magnetic 
surface  beyond  which  a  stochastic  region  appears  (cf .  Chapter  III) . 
Outside  that  magnetic  surface  the  particles  are  said  to  be  lost.  For 
simplicity  we  may  start  tracing  field  lines  with  at  initial  points  z  =-  0 
and  6-0  with  various  choice  of  r.  Suppose  the  magnetic  axis  is  at  r  - 
r  and  the  initial  point  for  the  outermost  magnetic  surface  is  at  r  - 
r^  .    Then  divide   the  interval  [r  ,  r.  ]  into  n  equal  parts  with  points 

r  ,  r,  ,  .  .  .  ,  r  -  r,  . 

o    1      '   n    L 
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Let   S  ,   S,  ,   . . . ,  S   be  the  magnetic  surfaces  traced  out  from  the 
o    J.         n 

initial  points  r  ,  r, r  .   Then  S   is  the  magnetic  axis  and  S   is 

o   i        n        o         ^  n 

the  boundary  surface.  Since  the  radial  variable  s  is  an  increasing 
function  of  the  magnetic  surfaces  and  has  the  value  0  at  the  magnetic 
axis  and  the  value  1  at  the  outermost  surface,  we  put  s(S,  )  —  (k/n)   for 

k  -  0,1 n,  where  a  is  a  positive  number.   We  see  that  s(S  )  -  0  and 

s(S^)  -  1. 

Suppose  the  position  of  a  particle  is  at  the  point  P  -  (r,fl,z)  in 
cylindrical  coordinates.  Then  the  magnetic  surface  S-  containing  P  lies 
between  two  magnetic  surfaces  S.  and  S.  .  or  outside  S  .  In  line  tracing 
of  the  magnetic  surfaces,  we  record  the  points  of  intersection  of  the 
field  line  and  the  cross  section  plane  v  -  v  .  If  we  follow  the  field 
line  sufficiently  many  times  around  the  torus,  an  invariant  closed  curve 
is  formed  on  the  plane  v  -  v  .  We  use  a  spline  fitting  routine  to 
approximate  the  closed  curve  and  interpolate  to  find  the  distance  from  a 
point  on  the  corresponding  surface  to  the  magnetic  axis  S  as  a  function 
of  its  poloidal  angle.  By  comparing  the  distance  of  the  point  P  to  the 
magnetic  axis  S  with  the  distance  functions  for  the  magnetic  surfaces, 
we  can  determine  between  which  two  magnetic  surfaces  S.  and  S.  ,  the 
point  P  lies,  or  whether  it  is  outside  the  outermost  magnetic  surface 
S  .  For  the  latter  case,  we  just  assign  a  value  greater  than  1  as  the 
radial  value  for  P.  For  the  former  case,  the  s  value  for  the  magnetic 
surface  S-  containing  P  is  in  the  interval  s(S.)  <  s(Sp)  <  s(S „  , ) .  We 
therefore  calculate  the  distance  of  P  to  the  surfaces  S.  and  S.  ,  in  the 
radial  direction  .  Using  interpolation  we  get  an  approximate  value  of  s 
for  Sp  and  hence  for  P. 
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Although  the  values  of  s  are  not  accurate,  our  method  seems  to 
guarantee  that  the  error  is  less  than  10%  if  n  -  10.  This  is  adequate 
to  estimate  the  confinement  time,  since  the  asymptotic  behavior  of  the 
Monte  Carlo  code  for  many  particles  averages  out  the  error  and  gives  at 
least  1  digit  of  accuracy. 

5 .   Scaling 

Since  the  drift  kinetic  equation  is  in  dlmensionless  form,  we  have 
to  specify  the  input  parameters  in  dimensionless  form  in  order  to  solve 
the  equation.  The  code  is  written  for  both  the  ion  and  electron  cases. 
Thus  when  we  pass  from  an  ion  run  to  a  corresponding  electron  run,  we 
have  to  adjust  the  parameters  in  a  prescribed  fashion. 

The  principal  input  parameters  for  the  Monte  Carlo  code  are  the 
gyroradius  of  the  particle  p^  and  the  collision  frequency  u .  The 
gyroradius  p.    in  units  of  the  plasma  radius  a  is  given  by 


.,  .^,  _Mv   ./2MT 

^^•^^)  ^L  -   qBa  "   qBa 


For  ions 


where  T  is  measured  in  kiloelectron  volts  and  B  in  tesla  and  a 
centimeters.  Typical  experimental  parameters  are  T  -  1  keV,  B  -  2  tesla 
and  a  -  20cm,  so  that  p.  =-0.01. 
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From  (5.15)  we  see  that  for  equal  ion  and  electron  temperatures  the 
gyroradlus  of  the  electron  will  be  smaller  than  that  of  the  ion  by  a 
factor  equal  to  the  square  root  of  the  ratio  of  the  mass  of  the  ion  to 
the  mass  of  the  electron.  Therefore  with  the  above  parameters,  the 
gyroradius  of  an  electron  is  roughly  0.0002. 

The  ion  collision  frequency  in  units  of  the  corresponding  cyclotron 
frequency  Is 


where   B   and  T  have   the   same   units   as  before  and  the  density  n  is 

3 
measured   in  particles   per   cm  .    For   convenience,   we   express   the 

collision   frequency   in  units   of   a  reference  collision  frequency  u 

When  T  -  10  keV  ,  B  =■  5T  and  n  =  2  x  lO''"^,  then  i/^^  -  0.6  x  lO"  . 


The  actual  collision  frequency  is  inversely  proportional  to  the 
square  root  of  the  mass  of  the  particle  [5].  Therefore  in  units  of  the 
corresponding  cyclotron  frequency,  which  is  inversely  proportional  to 
the  mass,  the  collision  frequency  becomes  proportional  to  the  square 
root  of  the  mass  of  the  particle.  Therefore  when  we  go  from  an  ion  run 
to  a  corresponding  run  for  an  electron,  w  is  reduced  by  a  factor  equal 
to  the  square  root  of  the  ratio  of  the  mass  of  the  ion  to  the  mass  of 
the  electron. 

The  principal  output  of  our  transport  code  is  the  confinement  time 
r,   which  is   expressed  in  units  of  the  reciprocal  of  the  cyclotron 
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frequency.  Thus  for  an  electron  case  the  time  is  reduced  by  a  factor 
equal  to  the  ratio  of  the  mass  of  the  ion  to  the  mass  of  the  electron. 

The  approximation  for  the  collision  operator  requires  that  i/6t  be 
small.  We  find  that  i/5t  -  0.01  is  adequate  for  the  ion  case.  For 
electrons  the  value  should  be  smaller,  say  0.005.  We  integrate  the 
drift  kinetic  equation  numerically  over  an  inteirval  [t  ,  t,],  where 

^d 

t,  -  t  -  7= 
1   o   yp^i/ 

Here  t .  is  1  for  the  ion  case  and  should  be  larger  for  the  electron 
case.  It  is  necessary  to  run  as  long  as  possible  to  get  adequate 
results.  In  order  to  provide  accurate  statistics  a  large  number  of 
particles  is  used.  To  optimize  the  computing  resources  and  the  result, 
we  use  128  particles,  which  should  give  a  reasonable  answer. 
Initialization  for  the  code  is  uniform  in  space  and  Maxwellian  in  the 
energy  profile. 

We  have  run  the  code  both  with  and  without  the  electric  potential 
in  both  the  ion  and  electron  cases.  The  results  suggest  that  the 
electric  potential  may  cause  island  overlapping  and  thus  account  for 
anomalous  transport  of  the  electrons.  This  confirms  the  contentions  we 
have  made  in  the  last  chapter. 


76 


VI.   THE  COMPUTER  PROGRAM 


We  have  combined  field  line  tracing,  evaluation  of  the  drift 
surface  island  width,  and  the  Monte  Carlo  method  of  calculating 
neoclassical  plasma  transport  in  one  code  call  FLDSMC .  Each  of  the 
options  is  based  on  integration  of  the  guiding  center  equations  (4.4). 
In  FLDSMC,  the  user  can  choose  to  run  one  of  the  three  options  by 
changing  the  control  parameter  NST.  For  field  line  tracing  NST  =■  1;  for 
the  drift  surface  island  width  NST  -  2;  and  for  the  Monte  Carlo 
calculation  NST  -  3.  In  Sections  1,  2  and  3  we  describe  some  of  the 
features  of  each  option  and  explain  how  to  operate  them. 

1.   Field  line  tracing 

The  first  option  allows  us  to  follow  magnetic  field  lines  or 
guiding  center  orbits  of  circulating  particles  around  the  torus  in  order 
to  study  the  structure  of  the  magnetic  surfaces  or  the  drift  surfaces. 
To  follow  field  lines  the  gyroradius  of  the  particle,  RADL,  is  set  to 
zero.   In  this  section  we  only  describe  magnetic  field  line  tracing. 

After  a  successful  run  tracing  the  magnetic  field  lines,  the 
magnetic  flux  surfaces  are  plotted  at  four  cross  sections  v  =  0.0,  0.25, 
0.50  and  0.75,  where  v  =  za/2n  is  the  normalized  toroidal  variable  with 
period  1.  The  points  of  the  graph  of  each  cross  section  are  stored  in 
separate  files.  The  user  is  allowed  to  enlarge  part  of  the  graph  at  any 
one  of  the  four  cross  sections  without  integrating  the  field  lines 
again.  A  circle  with  radius  equal  to  the  electron  gyroradius  is  plotted 
on   the  same  graph  for  comparison.   In  the  study  of  the  flux  surfaces  we 
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may  follow  the  field  lines  around  the  torus  at  least  ten  thousand  times 
in  order  to  have  an  adequate  picture  of  the  surface  on  the  scale  of  the 
electron  gyroradius  (cf.  Figs.  11-14).  The  user  can  continue  to  follow 
any  of  the  field  lines  by  setting  IFLAG  -  1. 

The  code  FLDSMC  uses  the  standard  ODE  solver  D02CHF  in  the  NAG 
FORTRAN  Library.  The  ODE  solver  D02CHF  uses  a  variable-order  variable- 
step  Adam's  method  to  integrate  a  system  of  first  order  ordinary 
differential  equations  until  a  user-specified  function  GCN  of  the 
solution  becomes  zero.  With  this  feature  we  can  stop  the  integration  of 
the  ODE  when  the  solution  for  the  variable  v  is  at  each  of  the  four 
cross  sections.  We  then  reset  the  function  GCN  and  restart  the 
integration  until  the  solution  for  v  is  at  the  next  cross  section.  The 
process  is  repeated  until  a  specified  number  of  field  periods  around  the 
torus  are  traversed.  This  feature  of  the  ODE  solver  provides  an 
accurate  solution  at  the  specified  cross  section  without  using 
interpolation  [18].  In  addition,  the  number  of  times  to  restart  the  ODE 
solver  is  reduced.  A  lot  of  computer  time  can  be  saved,  since 
restarting  the  ODE  solver  is  usually  expensive. 

The  magnetic  axis  in  the  general  3d  case  is  found  by  the  iteration 
method  described  in  Chapter  3.  In  the  2d  case  the  separatrix  is  also 
found  by  solving  the  system  of  equations  (3.3).  The  rotational 
transform  of  each  field  line  is  calculated  by  dividing  the  sum  of  the 
angles  through  which  the  field  line  rotates  about  the  magnetic  axis  each 
time  around  the  torus  by  the  number  of  field  periods  traversed. 

The  operation  of  the  code  is  relatively  simple.  The  user  must 
supply  the  number  of  field  lines  to  be  traced,  NORB,  the  number  of  field 
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periods  to  be  traversed  around  the  torus  for  each  field  line,  NVMAX,  and 
the  initial  position  of  each  field  line.  The  number  NORB  is  at  most  10 
and  can  be  set  by  modifying  the  read  statement  of  the  initial  position 
of  the  field  lines  in  the  main  routine.  The  parameter  LOG  controls  the 
enlargement  of  the  graph.  If  one  wishes  to  enlarge  part  of  the  graph, 
one  should  set  LOG  =  1  and  specify  the  coordinates  of  the  region  to  be 
enlarged  by  means  of  the  parameters  XMIN,  XMAX,  YMIN,  YMAX.  If  LOG  =-  0, 
no  enlargement  is  made . 

2 .   Drift  surface  island  width 

The  second  option  of  the  code  is  to  evaluate  the  island  width  of  a 
rational  drift  surface  at  the  v  =  0  cross  section.  The  gyroradius  of  the 
particle,  RADL,  the  initial  position  of  the  particle,  and  the  velocity 
pitch  angle  are  specified.  Typically  for  an  electron  RADL  =-  -0.0002 
and  for  an  ion  RADL  =  0.01.  The  minus  sign  indicates  an  electron  run. 
The  initial  position  is  required  to  be  near  to  the  desired  rational 
drift  surface.  An  approximate  initial  position  can  be  found  by  several 
trial  runs.  Once  an  approximate  initial  position  is  specified  the  code 
uses  an  iteration  method  to  find  the  elliptic  fixed  point  and  the  width 
of  the  corresponding  island  (cf.  Section  III.  3).  If  we  are  interested 
in  a  rational  drift  surface  with  the  rotational  transform  t  =-  n/ra,  the 
points  of  the  solution  on  the  specified  cross  section  are  saved  every  m 
times  around  the  torus.  This  can  be  accomplished  by  means  of  the  user- 
specified  function  GCN  in  the  ODE  solver.   The  electric  potential 

$  =  E0*cos(2^-kz) 
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in  cylindrical  coordinates  is  specified  by  the  parameters  EO ,  2  and  k.In 
order  to  monitor  the  calculation,  graphs  for  the  fixed  point  and  the 
magnetic  surfaces  of  the  island  are  plotted.  The  output  includes  the 
denominator  of  the  rotational  transform,  m,  the  position  of  the  elliptic 
fixed  point,  and  the  corresponding  island  width. 

3 .   The  Monte  Carlo  calculation 

The  last  option  of  the  code  calculates  the  neoclassical  confinement 
time  for  ions  and  electrons.  This  option  is  an  adaptation  of  the 
BETAMOC  code.  The  basic  difference  is  that  in  our  case  the  magnetic 
field  has  a  simple  analytical  representation  and  is  evaluated  internally 
in  the  code  so  that  an  equilibrium  solution  is  not  required.  However, 
we  have  to  pay  a  penalty  because  the  convenient  radial  variable  s  is  not 
available.  In  order  to  facilitate  the  run,  magnetic  surfaces  are 
calculated  by  field  line  tracing  at  60  values  of  v  before  the  run  and 
are  read  in  as  input.  The  number  of  particle  we  use  is  128.  If  more 
particles  are  needed,  we  have  to  change  the  number  NOL  in  the  PARAMETER 
statement  in  the  CLICHE  statements.  For  completeness  we  define  the 
input  parameters  in  a  sample  input  file.  More  details  can  be  found  in 
the  documentation  for  the  code  BETAMOC  [4]. 
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4.   Glossary  of  input  parameters 


CARDS  1,3 21 

CARD  2 
NS 

NPSI 

NPHI 

NMAG 

NEKIN 

CARD  4 
FREQ 

FREQE 

DELT 

TEND 
SEED 

NCOL 

CARD  6 
RADL 

EO 


NTEl,  NTE2 
TEl 

ACC 


FORTRAN  names   of   the  parameters.   Each  parameter 
occupies  8  columns. 


The   number   of   initial   values   in   the   radial 
direction.   FORMAT  18. 

The   number   of   initial  values   in   the   toroidal 
direction.   FORMAT  18. 

The   number   of   initial   values   in   the   poloidal 
direction.   FORMAT  18. 

The   number   of   initial   values   for   the  magnetic 
moment,    which   is   input   through   the   parameter 

n  =  V  /v.   FORMAT  18. 

11 

The   number   of   initial  values   for   the  kinetic 
energy.   FORMAT  18. 


The  collision  frequency  in  units  of  the  ion 
cyclotron  frequency  multiplied  by  0.6E-06.  FORMAT 
F8.4. 

The  collision  frequency  for  energy  scattering,  in 
the  same  units  as  FREQ.   FORMAT  F8 . 4 . 

The  time  interval  for  the  collision  operator. 
FORMAT  F8.4. 

The  time  of  the  run.   FORMAT  F8.4. 

The  initial  value  for  the  pseudorandom  number 
generator.   FORMAT  18. 

The  number  of  collisions  during  one  interval  of  the 
ordinary  differential  equation  solver.   FORMAT  18. 


The  gyroradius  in  units  of  the  plasma  radius.  RADL 
should  be  negative  for  electrons.   FORMAT  F8 . 4 . 

The  magnitude  of  the  electric  potential.  FORMAT 
F8.4. 

Parameters  in  the  background  temperature  profile. 
FORMAT  218.  F8 . 4  . 

The  exponent  in  the  accuracy  estimate  for  the 
ordinary  differential  equation  solver.  FORMAT 
F8.4. 
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CARD  8 
IC 
TLIM 

NST 

I  FLAG 

TITLE 

CARD  10 

ASPECT 

NU 


NV 

CARD  12 

BO.DEL-l.DELlO, 
DEL11.DEL21.DEL31 

CARD  14 

DEL12,DEL13,DEL22, 
DEL23,DEL32,DEL33 

CARD  16 

NT , NPT 


A  counter  used  to  monitor  printout.   FORMAT  18. 

The  estimated  total  CPU  time  for  a  run.  FORMAT 
F8 . 4 . 

The  control  parameter  of  the  run.  For  line  tracing 
NST  -  1;  for  drift  surface  island  width  NST  -  2; 
for  a  Monte  Carlo  calculation  NST  -  3.   FORMAT  18. 

An  indicator  to  continue  the  run.  For  a  normal  run 
IFLAG  -  0  and  to  continue  a  run  IFLAG  -  1.  FORMAT 
18. 

An  identifying  label  for  the  runs.   FORMAT  8A8 . 


The   ratio  of  the  major  radius  of  the  plasma  to  the 
minor  radius.   FORMAT  F8.4. 

The  number  of  points  used  in  the  poloidal  direction 
to  define  magnetic  surfaces.   FORMAT  18. 

The  number  of  points  used  in  the  toroidal  direction 
to  define  magnetic  surfaces.   FORMAT  18. 


The  coefficients  of  the  magnetic  field. 
FORMAT  6F8.4. 


More  coefficients  of  the  magnetic  field. 
FORMAT  6F8.4. 


These  numbers  define  the  maximum  number  of  terms 
used  in  each  of  the  double  sums  in  equation  (3.2). 
FORMAT  218. 


CARD 

18 

M,N 

RO 

ERR 

DRl, 

,DR2 

ET 


Parameters   used  to  define  the  form  of  the  electric 
potential.   FORMAT  218. 

Initial  guess  of  the  fixed  point.   FORMAT  F8.4. 

Tolerance  in  the  island  width.   FORMAT  F8.4. 

Step  sizes  in  calculating  the  island  width.   FORMAT 
2F8.4. 

Pitch  angle  of  the  velocity.   FORMAT  F8.4. 
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CARD  20 
NORB 
NSECT 
I  SECT 
LOG 

CARD  22 

XMIN.XMAX, 
YMIN.YMAX 

CARD  24 33 

STVl , STV2 , STV3 , 
NVMAX 


The  number  of  orbits  to  be  traced.   FORMAT  18. 

Number  of  sections  to  be  stored.   FORMAT  18. 

Section  to  be  enlarged.   FORMAT  18. 

A  control  parameter   for   enlarging  part   of  the 
graph.   To  enlarge  the  graph  LOG  -  1.   FORMAT  18. 


Coordinates   of   the  region  to  be  enlarged.   FORMAT 
4F8.4. 


The   initial  coordinates   for  line  tracing  and  the 
number  of  periods  to  be  traced.   FORMAT  3F8.4,  18. 
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5 .   Sample  run 


INPUT 


NS    NFS  I    NPHI    NMAG     NEK 
4       2       4       2       2 


FREQ    FREQE     DELT    TEND     SEED    NCOL 
4.0000  4.0000   0.0025   0.0001   507.00       1 


RADL  EO  NTEl  NTE2  TEl     ACC 

0.01000  0.0100  0.0000  0.0000  0.9900   5.0000 

IC  TLIM  NST  I  FLAG  TITLE 

1  7200.  1  0  LINE  TRACING 


ASPECT      NU      NV 
2.0000      60      60 


DELO    DEL-1    DELIO    DELll    DEL21    DEL31 
1.0000   0.1000   0.0000   0.1850  -0.1100   0.0000 


DEL12    DEL22    DEL32   DEL13   DEL23   DEL33 
0.0000   0.0000   0.0000   0.0000   0.0000   0.0000 


NT     NPT 
2       1 


M       N      RO      ERR     DRl      DR2       ET 
17      10   0.8000  0.00001   0.0001   0.0005   0.8000 


RESULT  OF  FIELD  LINE  TRACING 

ORBIT     INITIAL  POSITION       IOTA     PERIODS 
1      0.8000  0.0  0.0     0.2795     100.0 

COMPUTING  TIME  =       5. 7  3 
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6 .   Listing  of  the  code 

TABLE  OF  SUBROUTINES 
C 
C 
C      0.  CLICHE  -  COMMON  BLOCKS  AND  DIMENSION  PARAMETERS        87 

C 

C       1.  FLDSMC  -  MAIN  CONTROL  PROGRAM  88 

C 

C       2.  FLT     -  FIELD  LINE  TRACING  91 

C 

C       3.  LOR     -  CONTROL  OF  THE  PLOTTING  FOR  LINE  TRACING       92 

C 

C       4.  ORBPLT  -  MONITORING  OF  THE  LINE  92 

C 

C       5.  ROTNU   -  EVALUATION  OF  THE  ROTATIONAL  TRANSFORM        93 

C 

C      6.  CROSEC  -  FINDING  THE  PLANE  OF  INTERSECTION  94 

C 

C       7.  PPLOT   -  PLOTTING  LINES  IN  DIFFERENT  CROSS -SECTIONS     95 

C 

C       8.  LARGEP  -  ENLARGEMENT  OF  ONE  SECTION  95 

C 

C       9.  LOCAL   -  ENLARGING  PART  OF  THE  GRAPH  96 

C 

C      10.  CIRCLE  -  DRAWING  A  CIRCLE  97 

C 

C     11.  MAXIS   -  DETERMINATION  OF  THE  MAGNETIC  AXIS  IN  2D      97 

C 

C     12.  PSIP    -  STREAM  FUNCTION  97 

C 

C      13.  FMAXIS  -  DETERMINATION  OF  THE  MAGNETIC  AXIS  IN  3D       98 

C 

C      14.  DISLD   -  DRIFT  SURFACE  ISLAND  WIDTH  99 

C 

C      15.  FIXPT   -  FINDING  THE  ELLIPTIC  FIXED  POINT  101 

C 

C      16.  WIDTH   -  WIDTH  OF  AN  ISLAND  102 

C 

C      17.  GRAPH   -  GRAPH  FOR  THE  ISLAND  104 

C 

C      18.  ORBIT   -  GUIDING  CENTER  DIFFERENTIAL  EQUATIONS         104 

C 

C      19.  EVAL    -  EXPONENTIAL  DECAY  FUNCTIONALS  107 

C 

C     20.  COLL    -  COLLISION  OPERATOR  109 

C 

C      21.  BFIELD  -  EVALUATION  OF  MAGNETIC  FIELD  109 

C 

C      22.  GCN2    -  SPECIFIED  FUNCTION  FOR  ORBITS  110 

C 

C      23.  GCN     -  SPECIFIED  FUNCTION  :--0R  TRANSPORT  111 

C 

C      24.  FCN     -  RIGHT-HAND  SIDE  OF  ODE  HI 

C 

C      25.  BESSEL  -  EVALUATION  OF  MODIFIED  BESSEL  FUNCTIONS       112 

C 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


26.  INIT 

27.  ENIN 

28.  TEMP 

29.  SURF 

30.  SFCT 

31.  SSS 

32.  TRANSF  - 

33.  SPLIF   - 

34.  INTPL   - 

35.  TPLOT   - 

36.  AVER 

37.  DENSITY- 

38.  SLOPE   - 

39.  PLOTB   - 

40.  ELPOT   - 


RANSET 
RANF 
SECOND 
D02CHF 


KEEP80 

FR80ID 

MAP 

FRAME 

GAXIS 

SETLCH 

CRTBCD 

SETCRT 

VECTOR 

POINT 


INITIALIZATION  OF  VARIABLES 

ENERGY  INITIALIZATION 

TEMPERATURE  PROFILE 

DEFINE  MAGNETIC  SURFACES 

EVALUATION  OF  THE  RADIAL  VARIABLE 

RADIAL  FUNCTION 

TRANSFORMATION  INTO  CYLINDRICAL  COORDINATES 

GENERAL  PURPOSE  CUBIC  SPLINE 

GENERAL  PURPOSE  INTERPOLATION 

PLOT  OF  LOSS  RATE  AND  ORBITS 

PARTICLE  FUNCTION  AVERAGE 

PARTICLE  DISTRIBUTION 

ESTIMATION  OF  LOSS  RATE 

DRIVER  FOR  PLOT  ROUTINES 

EVALUATION  OF  THE  ELECTRIC  POTENTIAL 

INDEX  OF  LIBRARY  ROUTINES 


112 
113 
114 
114 
116 
117 
118 
118 
119 
120 
125 
126 
126 
127 
128 


ESTABLISH  RANDOM  NUMBER  SEED 
OBTAIN  RANDOM  NUMBER 
CUMULATIVE  CPU  TIME 
INTEGRATE  A  SYSTEM  OF  FIRST  ORDER 
ORDINARY  DIFFERENTIAL  EQUATIONS  WITH 
INITIAL  CONDITIONS  USING  A  VARIABLE- 
ORDER  VARIABLE -STEP  ADAMS  METHOD 
D02CHF  IS  IN  THE  NAG  LIBRARY 
SAVE  PLOT  FILE  FOR  LATER  USE 
MAKE  AN  FR80  FILE 
TRANSLATE  AND  SCALE  PLOT 
ADVANCE  FRAME  95 

DRAW  AND  LABEL  AXES 
PLOT  IN  USER'S  COORDINATES 
WRITE  LABELS  ON  AXES 
SET  INITIAL  POSITION  FOR  PLOT 
DRAW  A  LINE 
DOT  A  POINT 


90 

90,109 

90 

92,100,102 


92 
92 

95,104,122 

96,116,121 

95,96,104 

121,127 

121,127 

121,127 

128 

95,97 


THE  PRECOMP  UTILITY  PROCESSES  PARAMETERS,  CLICHES,  AND  ENDCLICHES 
THE  OUTPUT  OF  PRECOMP  CAN  BE  SUBMITTED  TO  THE  FORTRAN  COMPILER. 
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CLICHE  NAMEB 

PARAMETER  (NOL=128,  N0L4=4*N0L) 

ENDCLICHE 

CLICHE  NAMEl 

USE  NAMEB 

COMMON  /BF/DBB(NOL,3,3) , BB(NOL, 3) ,BVAL(NOL) 

ENDCLICHE 

CLICHE  NAME2 

COMMON  /BP/NT,NPT,DEL(4,3) , BO , DELIO , DELMl 

ENDCLICHE 

CLICHE  NAME 3 

USE  NAMEB 

COMMON  /MC/  TJ(NOL),  SJ(NOL),  FREQ ,  RFREQ ,  SFREQ,  SDELT,  TSE, 
1  FTAU,  Tl.  T2,  NLOST,  NCOL,  NST,  ICOL,  TSEl,  NPRT 

ENDCLICHE 

CLICHE  NAME4 

USE  NAMEB 

COMMON/INI/  PHIO(NOL) ,PSI0(NOL) ,S0(NOL) ,ETA0(NOL) ,AMAG0(NOL) , 
1  NPAR, EKING (NOL),R0(NOL) ,U0(NOL) ,TE1 ,NTE1 .NTE2 ,T0 

ENDCLICHE 

CLICHE  NAME5 

COMMON  /CPL/  ALP,  IC,  10,  TLIM,  ACC ,  NTER,  KTER, 
1  NOR,  TIMF 

ENDCLICHE 

CLICHE  NAME6 

COMMON  /POL/  RADL,  EC,  DELT,  TEND,  EG,  XEO ,  EGC,  PCO, 
1  TSOLl,  El,  E2,  E3,  lEL 

ENDCLICHE 

CLICHE  NAME? 

DIMENSION  XSURF(51G) ,YSURF(51G) 

ENDCLICHE 

CLICHE  NAMES 

COMMON  /SURFP/SURFR(51G,62) , SURFZ(51G , 62) ,RA(62) ,ZA(62) , 
1  NU,NV,HU,HV,NI 

ENDCLICHE 

CLICHE  NAME9 

COMMON/EVAL/  FAV , GAV , ETAV , ETDV , PAV , SAV , SDV , DENER , NLEFT , 
1  FACS,TEAV,TEMAX,TEMIN 

ENDCLICHE 

CLICHE  NAMEIO 

USE  NAMEB 

COMMON  /SOL/  SS(NOL) ,PSIT(NOL) ,PHIT(NOL) ,ETA(NOL) ,AMAG(NOL) , 

1  ENER(NOL) ,ETAMA(NOL) ,ETAMI(NOL) .PESl(NOL) ,PES2(N0L) ,PES3(N0L) , 

2  WORK(NOL4,18) ,EKIN(NOL) ,EKMAX(N0L) ,EKMIN(N0L) ,PES4(N0L) , 

3  PES5(N0L) ,PES6(N0L) ,P0T(N0L) ,POTR(NOL) ,POTU(NOL) ,P0TZ(N0L) 
ENDCLICHE 

CLICHE  NAMEll 

USE  NAMEB 

DIMENSION  Y(N0L4) ,YPRIME(N0L4) 

ENDCLICHE 

CLICHE  NAME12 

COMMON/PLOT/  PL(700G , IG) , IN(5) 

ENDCLICHE 

CLICHE  NAME13 

COMMON  /PIS/  PI2,  PI 

ENDCLICHE 

CLICHE  NAME14 

DIMENSION  UR(NOL) ,SR(NOL) 


87 


ENDCLICHE 

CLICHE  NAME15 

COMMON  /CONT/  YC(N0L4) 

ENDCLICHE 

CLICHE  NAME16 

COMMON  /WDPAR/  M.N.RO, ERR, DRl , DR2 

ENDCLICHE 

CLICHE  NAME17 

COMMON  /CONST2/VPL(4) .R0T(4) ,NPER(4) , IPL(4)  ,  IPL0T(4) 

COMMON  /C0NST3/RDMA(4) ,RDMI (4) , IC0UNT(4 , 20) ,ROTF(4,20) 

COMMON  /CONST4/LOST(20) ,NSTEP(20) 

COMMON  /CONST5/FPP(20) ,TIME 

COMMON  /CONST7/RAX,RX(4) ,RY(4) ,DX1,DY1,T0L 

COMMON  /VAR4/NVMAX(20) ,STV(3,20) 

COMMON  /VAR5/N0RB,NSECT,IFLAG 

COMMON  /VARII/DV,FINV 

COMMON  /IOTA/  ROTN(20) 

ENDCLICHE 


PROGRAM  FLDSMC 

THIS  MAIN  ROUTINE  CALLS  OTHER  SUBROUTINES  AND  CONTROLS  INPUT-OUTPUT 

USE  NAME2 

USE  NAME3 

USE  NAME4 

USE  NAMES 

USE  NAME6 

USE  NAMES 

USE  NAME9 

USE  NAMEll 

USE  NAME12 

USE  NAME13 

USE  NAME15 

USE  NAME16 

USE  NAME17 

DIMENSION  NW(20)  ,TITLE(4) 

DATA  I PLOT/0, 0,0,0/ 

CALL  LINK  ("UNIT25-(INDATA, OPEN) ,UNIT6-(0UTPUT, CREATE, TEXT) , 

1  UNIT002-(TAPE2, CREATE,     ) ,UNIT003-(TAPE3 , CREATE.       ) 

2  UNIT007-(TAPE7, CREATE,     ) ,UNIT001-(TAPE1 , CREATE,     ), 

3  UNIT004-(TAPE4, CREATE,     ),READ25, 

4  PRINT6//") 
READ  (25.600) 

READ  (25,610)  NS ,NPSI ,NPHI ,NMAG,NEKIN 

READ  (25.600) 

READ  (25,620)  FREQ, FREQE, DELT, TEND, SEED, NCOL 

READ  (25,600) 

READ  (25,630)  RADL, EO ,NTE1 ,NTE2 ,TE1 ,ACC 

READ  (25.600) 

READ  (25.640)  IC ,TLIM,NST, IFLAG, TITLE 

EC-0 . 5*RADL*RADL 

PRINT  900 

PRINT  910.NS,NPSI,NPHI,NMAG,NEKIN 

PRINT  920 . FREQ , FREQE , DELT , TEND , SEED , NCOL 

IEL-1 

IF(RADL.LE.O.O)  IEL--1 
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10 


20 


PRINT  930 , RADL , EO , NTEl , NTE2 , TEl , ACC 

RADL-ABS(RADL) 

PRINT  940, IC.TLIM.NST.IGLAG, TITLE 

RFREQ-0.6E-06 

S FREQ-FREQ*RFREQ 

EFREQ-FREQE*RFREQ 

SDELT-DELT/AMAXl ( SFREQ , EFREQ) 

TSE-TEND*SQRT(FREQ)/(10.*SFREQ*SQRT(RADL)) 

FTAU-SDELT*SFREQ 

FTAUE-SDELT*EFREQ 

SDELT-SDELT*NCOL 

DEFINE  CONSTANTS 

PI-3. 1415926535898 

PI2-2.0*PI 

READ  IN  PARAMETERS 

READ  (25,600) 

READ  (25,650) 

ALP-l.O/ASP 

READ  (25,600) 

READ  (25,660) 

READ  (25,600) 

READ  (25,660) 


IN  MAGNETIC  FIELD 
ASP,NU,NV 

BO,DELM1,DEL10,DEL(1,1),DEL(2,1) ,DEL(3,1) 

2).DEL(1,3) 


READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

READ 

DO  10 

READ 


(25,600) 
(25,670) 
(25,600) 
(25,680) 
(25,600) 
(25,690) 
(25,600) 
(25,710) 
(25,600) 
1-1,10 
(25,700) 


DEL(1,2),DEL(2,2),DEL(3 
,DEL(2,3) ,DEL(3,3) 

NT.NPT 

M,N,R0,ERJI,DR1,DR2,ET 

NORB,NSECT,ISECT,LOC 

XMIN , XMAX , YMIN , YMAX 


STV(1,I) ,STV(2,I) ,STV(3,I) ,NVMAX(I) 


CONTINUE 
PRINT  950 
PRINT 
PRINT 
PRINT 


960, 
970, 
980, 


ASP,NU,NV 

BO,DELM1,DEL10,DEL(1,1) ,DEL(2,1) ,DEL(3,1) 

DEL(1,2) ,DEL(2,2) ,DEL(3,2) ,DEL(1,3) 
,DEL(2,3) ,DEL(3,3) 

NT,NPT 

M,N,R0,ERR,DR1,DR2,ET 
20,  110,  120 

TRACING 


30 


40 


PRINT  990, 

PRINT  995, 

IF(NST-2) 

CALLL  LINE 

PRINT  800 

NLEFT-1 

TOL-10**(-ACC-4.0) 

FIND  THE  MAGNETIC  AXIS 

CALL  FMAXIS(IFAIL,R1) 

PRINT  820 

REWIND  1 

IF(IFL.\G.GE.O)  GO  TO  30 

GO  TO  100 

IF(IFLAG.EQ.O)  GO  TO  50 

REWIND  2 

DO  40  L-1,N0RB 

READ(2)  X1,X2,X3,STV(1,L)  .STV(2,L)  ,STV(3,L)  ,NW(L)  ,R0TN(L) 

NVMAX(L)-NVMAX(L)+NW(L) 
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READ(2)  (IPL0T(I),I=1,4) 
50   REWIND  2 

NOR-NORB 

DO   70  L-1,  NOR 

Y(l)    -STV(1,L) 

Y(2)    -STV(2,L) 

Y(3)    -STV(3,L) 

CALL  FLT(Y,LOST(L) ,NSTEP(L) ,NVMAX(L) , L) 

WRITE(2)  STV(1.L),STV(2,L),STV(3,L),Y(1),Y(2),Y(3),NVMAX(L) 
1         ROTN(L) 

NORB-L 

X-Y(2)/Y(3) 

PRINT  830,L,STV(1,L) ,STV(2.L) ,STV(3,L) .ROTN(L) , FPP(L) .X 

IF(LOST(L) .EQ.l)  GO  TO   80 
70   CONTINUE 
80   DO  90  1-1,6 
90   FPP(I)--1.0 

URITE(2)  (FPP(I) ,1-1,6) ,NVMAX(1) ,R0TN(1) 
100   CALL  LOR 

STOP  7  7 
C     DRIFT  SURFACE  ISLAND  WIDTH 
110   Rl-RO 

CALL  DISLD(R1,ET) 

STOP  88 
C     MONTE  CARLO  CALCULATION  OF  TRANSPORT 
120   ISEED-ABS(SEED) 

CALL  RANSET(ISEED) 

ICOL-0 

IF(NST.LE.O)  GO  TO  160 

REWIND  4 

DO  130  1-1,20000 

READ  (4)  (PL(I.J) ,J-1,10) 

IF(PL(I,1) .LT.0.0)  GO  TO  140 
130   CONTINUE 
140   READ  (4)  X1,X2,TSE1,NM1,NM2,IC0L,NM4,NM5.NM6 

TSE-TSE+TSEl 

DO  150  I-1,IC0L 

Xl-RANF(l) 

X2-RANF(1) 
150   CONTINUE 
160   CONTINUE 

GALL  ORBIT (NS , NPHI , NPS I , NMAG , NEKIN , FTAUE) 

PRINT  992 

CALL  TPLOT(YC) 

STOP  99 
600   FORMAT(10A8) 
610   FORMAT(5I8,) 
620   FORMAT(5F8.4,I8) 
630   F0RMAT(2F8.4,2I8,F8.4) 
640   FORMAT(I8,F8.4.I8,I8.4A8) 
650   FORMAT(F8.4,2I8) 
660   F0RHAT(6F8.4) 
670   F0RMAT(2I8) 
680   FORMAT(2I8,5F8.4) 
690   F0RMAT(4I8) 
700   FORMAT(3F12.8,ir2) 
710   FORMAT (4F8. 4) 
800   FORMAT (///9X, 'RESULT  OF  FIELD  LINE  TRACING') 
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810   FORMAT(//13X,'2D  MAGNETIC  AXIS  R  =  ',  F7 . 3// 

1        13X,'3D  MAGNETIC  AXIS  R  =  ',  F7.3) 
820   FORMAT(//9X,5HORBIT,4X, 'INITIAL  POSITION' ,4X,'   I0TA',4X, 

I        'PERIODS' ,4X, 'U/V  RATIO'//) 
830   FORMAT(I0X,I2,4X,F8.4,2F5.1,4X,F7.4.3X,F7.1,4X,F8.3) 
840   F0RMAT(//,6X, 'COMPUTING  TIME  =  ',F10.2/) 

900   F0RMAT(///36X,5HINPUT,/) 

910   FORMAT (//,11X,2HNS,4X,4HNPS I. 4X,4HNPHI,4X,4HNMAG,5X,3HNEK,4X, 

1  /5X,5I8) 
920   F0RMAT(// , 9X , 4HFREQ , 3X , 5HFREQE , 4X , 4HDELT , 4X , 4HTEND , 4X , 4HSEED , 

1  4X,4HNC0L,/  5X,4F8.4,F8.2,I8) 
9  30   FORMAT ( // , 9X , 4HRADL , 6X , 2HE0 , 4X , 4HNTE1 , 4X , 4HNTE2 , 5X , 3HTE1 , 5X , 

1  3HACC/5X.F8.5,4F8.4,F8.4) 
940   FORMAT (/,11X,2HIC,4X,4HTLIM,5X,3HNST,3X,5HI FLAG, 

1  3X, 5HTITLE/5X , 218 , F8 . 1 . 218 , 4A8) 
950   FORMAT (////27X, 'PARAMETERS  OF  THE  MAGNETIC  FIELD'/) 
960   F0RMAT(//11X, 'ASPECT' ,6X, 'NU' ,6X, 'NV'/9X, F8 .4, 218) 
970   F0RMAT(//13X, 'DELO' ,3X, 'DEL-1' ,3X, 'DELIO' ,3X, 'DELll' ,3X, 

1       'DEL21' ,3X, 'DEL31'/9X,6F8.4) 
980   F0RMAT(//12X, 'DEL12' ,3X, 'DEL22' ,3X, 'DEL32' ,3X, 'DEL13' ,3X, 

1       'DEL23' ,3X, 'DEL33'/9X,6F8.4) 
990   F0RMAT(//15X, 'NT' ,5X, ' NPT' /9X, 218) 

992   FORMAT (//////, 3 5X," OUTPUT",//) 

994  FORMAT (//I OX, 'COMPUTING  TIME  =  ',F8.2) 

995  FORMAT (//16X. 'M' , 7X, ' N' , 6X, 'RO' , 5X. 'ERR' ,5X. 'DRl' ,5X, 'DR2' , 
1   6X, 'ET'/9X,2I8,5F8.5) 

END 


SUBROUTINE  FLT ( Y , LOSTl , NSTEPl , NVM , L) 

PERFORM  FIELD  LINE  TRACING 

USE  NAME13 

USE  NAME17 

DIMENSION  Y( 3) ,  WORK(3,18) 

EXTERNAL  FCN,GCN2 

LOST 1-0 

H-20.0 

DV-1 . 0/FLOAT(NSECT) 

VMAX-NVM 

NVSTEP-NSECT*NVM 

NSTEPl-1 

Tl-0.0 

T2-0.0 

HMAX-0.0 

NODE-3 

IF(IFLAG.GT.O)  GO  TO  10 

FINV-DV 

DX1-Y(1)-RAX 

DY1=0.0 

WRITE(l)  Y(l) ,Y(2) ,Y(3) 

GO  TO  20 
10   FINV-Y(3)+DV 

DX1-Y(1)*C0S(PI2*Y(2))-RX(K) 

DY1-Y(1)*SIN(PI2*Y(2))-RY(K) 
20   T2-T2+H 

IFAIL-1 
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CALL  D02CHF(T1,T2, NODE, Y,T0L,1.HMAX,FCN,GCN2, WORK, IFAIL) 

T2-T1 

IF(IFAIL.NE.O)  GO  TO  30 

FINV-FINV+DV 

WRITE(l)  Y(1),Y(2),Y(3) 

NSTEPl-NSTEPl+1 
30   CONTINUE 

IF  (Y(l) .LT.O.O.OR.Y(l) .GT.5.0)  GO  TO  UO 

IF(NSTEPl.GE.NVSTEF)  GO  TO   50 

IF(Y(3) .GE.VMAX)  GO  TO   50 

GO  TO  20 
C     ORBIT  LOST 
40   LOSTI-1 
50   FPP(L)-Y(3) 

RETURN 

END 


SUBROUTINE  LOR 
:     INITIATE  PLOTTING 

USE  NAMEI3 

USE  NAME17 

DATA  VPL/0.0,0.25,0.50.0.75/ 

CALL  SECOND(TIMl) 

SCAL-2.0 

CALL  KEEP80(1,3) 

CALL  FR80ID 
5   DO  10  K-1,4 

KK-IO+K 

REWIND  KK 
10   CONTINUE 

REWIND  1 

PRINT  900 

CALL  ORBPLT 

IF(IFLAG.LE.O)  CALL  ROTNU 

CALL  PPLOT(SCAL) 

CALL  LARGER (l.SCAL) 

IF(LOC.GT.O)  CALL  LOCAL( ISECT.NP) 

CALL  SEC0ND(TIM2) 

TIME-TIM2-TIM1 

PRINT  910,  TIME 

RETURN 
900   FORMAT (//4X, 'RESULT'//) 
910   F0RMAT(//,6X, 'COMPUTING  TIME  -  '.F7.2/) 

END 


SUBROUTINE  ORBPLT 

MONITOR  THE  LINE  THAT  FALLS  INTO  4  CROSS -SECTIONS 

USE  NAME13 

USE  NAME17 

DATA  IPL/0, 0,0,0/ 

IF(IFLAG.EQ.O)  GO  TO  30 

DO  20  K-1,4 
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DO  10  I=1,IPL0T(K) 

READ(K+10,END-20,ERR-20)  A,B 
10   CONTINUE 
20   CONTINUE 
30   DO   70  L-1,N0RB 

DO  40  NL-1,NSTEP(L) 

READ(l,END-50,ERR-50)  R2,U2,V2 

CALL  CR0SEC(V2,K) 

UP-PI 2*U2 

ZP-R2*SIN(UP) 

RP-R2*C0S(UP) 

WRITE(10+K)  RP.ZP 

IPL(K)-IPL(K)+1 
40   CONTINUE 
50   CONTINUE 

DO  60  K-1,NSECT 

ICOUNT(K,L)-IPL(K) 
60   CONTINUE 
70   CONTINUE 

WRITE(2)  (IPL(I) ,1-1,4) 

RETURN 

END 


SUBROUTINE  ROTNU 

FIND  THE  ROTATIONAL  TRANSFORM  FOR  EACH  MAGNETIC  LINE 

USE  NAME13 

USE  NAME17 

DIMENSION  SUMTH(4) ,DIFR1(4) ,DIFR2(4) ,DIFU1(4) ,DIFU2(4) 

PRINT  900 

PRINT  910, VPL 

DO  10  K=1,NSECT 

KK=K+10 

REWIND  KK 
10   CONTINUE 

IF(IFLAG.EQ.O)  GO  TO  50 

DO  30  K-1,4 

DO  20  I-1,IPL0T(K) 

READ(K+10,END-30,ERR=30)  A,B 
20   CONTINUE 
30   CONTINUE 

DO  40  1-1,4 
40   IPLOT(I)-IPLOT(I)+IPL(I) 

WRITE(3)  (IPLOT(I) , 1-1,4) 
50   PRINT  920,  (RX(K)  ,RY(K),K=-1,NSECT) 

PRINT  930 

DO  110  L-1,N0RB 

DO  100  K-1,NSECT 

NPER(K)  -1 

SUMTH(K)=0.0 

KK=K+10 

M-1 

IF(L.GT.l)  M=IC0UNT(K,L-1)+1 

DO  80  I-M,ICOUNT(K,L) 

READ(KK,END=90,ERR=90)  RP,  ZP 

DIFR2(K)-RP-RX(K) 
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55 

60 
70 


80 
90 


100 

110 

900 
910 
920 
930 
9A0 


DIFU2(K)-ZP-RY(K) 
IF(NPER(K) .EQ.l)  GO  TO  70 
BB-DIFR2(K)**2+DIFU2(K)**2 
AA-DIFR1(K)**2+DIFU1(K)**2 
A-SQRT(BB) 
B-SQRT(AA) 

D0TPD-DIFR1(K)*DIFR2(K)+DIFU1(K)*DIFU2(K) 
CR0PD-DIFR1(K)*DIFU2(K)-DIFU1(K)*DIFR2(K) 
IF(A*B.EQ.O.O)  GO  TO  55 
X-DOTPD/(A*B) 

IF  (ABS(X) .GT.1.0)  GO  TO  55 
THETA-ACOS(X) 

IF(CROPD.LT.O.O)  THETA-PI2 -THETA 
GO  TO  60 
THETA-0 . 0 

IF(X.LT.O.O)  THETA-0.5*PI2 
SUMTH(K)-SUMTH(K)+THETA 
NPER(K)-NPER(K)+1 
DIFR1(K)-DIFR2(K) 
DIFU1(K)-DIFU2(K) 
CONTINUE 
CONTINUE 

IF  (NPER(K) .GT.l)  NPER(K)=NPER(K) - 1 
ROT(K)  -  SUMTH(K)/FL0AT(NPER(K))/PI2 
R0TF(K,L)-R0T(K) 
CONTINUE 

PRINT  940,L,VPL,ROT,NPER 
RETURN 

FORMAT ( 2 2X, 'REFERENCE  POINT  IN  EACH  CORSS- SECTION  PLANE'//) 
F0RMAT(9X, 'SECTION' ,4F15. 5) 

F0RKAT(9X, 'REF.  PT' , 4(2X , ' ( ' , F5 . 2 , ' , ' , F5 . 2 , ' ) ' ) ) 
F0RMAT(//29X, 'ROTATION  NUMBER  OF  EACH  LINE'//) 
F0RHAT(//9X, 'LINE  ' ,I3/9X, 'SECTION' ,4F15.5/ 
+  9X.'R0T#   '  ,4F15.5/9X, ';s  PTS   ',4115) 
END 


10 


900 


SUBROUTINE  CROSEC(V , KVP) 

FIND  THE  SECTION  THAT  A  POINT  FALLS  IN 

USE  NAME13 

USE  NAME17 

DIMENSION  VP(4) 

DATA  VP/O. 0,0. 25, 0.50, 0.75/ 

VMAX-AMOD(V,1.0) 

VMIN-AMOD(V,1.0) 

IF(VMIN.GT.0.9)  VMIN=VMIN- 1 . 0 

VMAX-VMAX+0.001 

VMIN-VMIN-0.001 

DO  10  K-1,NSECT 

KVP-K 

IF(VMAX.GE.VP(K) . AND . VMIN . LE. VP(K) )  RETURN 

CONTINUE 

PRINT  900 

RETURN 

FORMAT (//9X. 'NO  MATCH  OF  THE  PLANE') 

END 
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10 


20 
30 


910 
920 
900 


SUBROUTINE  PPL0T(SCAL2) 

TO  PLOT  THE  LINES  IN  4  DIFFERENT  CROSS- SECTIONS 
USE  NAME2 
USE  NAME13 
USE  NAME17 

DIMENSION  Xl(4) ,X2(4) ,Y1(4) ,Y2(4) 
DATA  Xl/0. 121, 0.458, 0.121. 0.458/ 
DATA  Yl/0. 586, 0.586, 0.253. 0.253/ 
FACl-0.269 
SCAL1--1.0*SCAL2 
DO  30  K-1,NSECT 
X2(K)-X1(K)+FAC1 
Y2(K)-Y1(K)+FAC1 

CALL  MAP(SCAL1,SCAL2,SCAL1,SCAL2,X1(K) ,X2(K) ,Y1(K) ,Y2(K)) 
CALL  LINE  (SCAL1,0.0,SCAL2, 0.0,0) 
CALL  LINE  (0.0,SCAL1,0.0.SCAL2.0) 
KK-IO+K 
REWIND  KK 

READ(KK,END-20,ERR-20)  RP , ZP 
CALL  P0INT(RP,ZP,1) 
GC  TO  10 
CONTINUE 
CONTINUE 

CALL  MAP(0.0,16.0,0.0,16.0,X1(1),X2(2).Y1(3),Y2(2)) 
CALL  SETLCH(-1.0,0.0,1,0,1,0) 
ASP=1.0/ALP 
WRITE(100,900) 

WRITE(100,910)  ASP,DELM1,DEL(1,1) 
WRITE(100,920)  (kOTN(I) . 1=1 ,N0RB) 
CALL  FRAME 
RETURN 
F0RMAT(/2X, 
FORMAT (/2X, ' 
FORMAT (//2X, 
1        ,  '  .50 
END 


ASPECT  =' ,F5.2,3X,7HDEL-1  =,F5 
IOTA  =' ,10(2X,F4.3)) 
'MAGNETIC  SURFACES  OF  A  HELIAC 
75') 


2,3X,7HDEL11  =,F5.2) 
FIELD  AT  V  =  .00, .25, ' 


10 


SUBROUTINE  LARGEP( ISECT , SCAL2) 

PLOT  THE  ISECT  SECTION  ON  A  SMALLER  SCALE 

USE  NAME13 

USE  NAME17 

DIMENSION  DIV(2) 

DATA  XI, X2,Yl,Y2/0. 12 1,0. 706, 0.253, 0.838/ 

SCAL1--1.0*SCAL2 

KK-10+ISECT 

REWIND  KK 

CALL  MAP ( SCALl , SCAL2 , S CALI , SCAL2 , XI , X2 , Yl , Y2 ) 

DIV(l)=-3 

DIV(2)-2. 

CALL  GAXI S(SCAL1,0.0,S  CAL2 , 0 . 0 , 0 , 1 , 0  .  "  1 2 "  ,  2  ,  DI V ) 

CALL  GAXIS (0 . 0 , SCALl ,0.0, SCAL2 , 0 , 1 , 1 , " 12 " , 2 , DIV) 

CONTINUE 

READ(KK,END=20,ERR=20)  RP,ZP 

CALL  POINT(RP,ZP, 1) 

GO  TO  10 
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20   CONTINUE 

CALL  SETLCH(SCAL1,SCAL1.1,0,1,0) 

WRITE(100,900)  RAX 
900   FORMAT (//IX, 'MAGNETIC  AXIS   R-  ',F5.3) 

WRITE(100,910)  (ROTN(I) ,I-1,N0RB) 
910   FORMAT (/ IX, 'IOTA  FOR  EACH  SURFACE  STARTING  FROM  MAGNETIC  AXIS:'// 
I       .10F6.3) 

CALL  FRAME 

RETURN 

END 


SUBROUTINE  LOCAL( ISECT , NP) 
C     ENLARGE  PART  OF  THE  GRAPH  IN  ISECT  SECTION  DEFINED  BY 
C      (XMIN.XMAX)  X  [YMIN.YMAX] 

USE  NAME13 

USE  NAME17 

DIMENSION  XA(2) ,YA(2) 

DATA  XI , X2 , Yl , Y2/0 .130,0.650,0.14,0.66/ 

XA(1)-XMIN 

XA(2)-XMAX 

YA(1)-YMIN 

YA(2)-YMAX 

RAD-0.0002 

PRINT  900,XA(1),XA(2),YA(1),YA(2) 

CALL  MAP(XA(1) ,XA(2) ,YA(1) ,YA(2) ,X1 ,X2 , Yl . Y2) 

CALL  GAXIS(XA(1) , YA( 1) ,XA(2) , YA(1) , 0 , 0 , 0 , "F8 . 4" , 3 .XA) 

CALL  GAXIS(XA(1) , YA(1) ,XA(1) , YA(2) , 0 , 0 , 1 , "F8 .4" , 3 , YA) 

KK-ISECT+10 

REWIND  KK 

XMAX--99.0 

DO  30  L-1,N0RB 

XMA--99.0 

XMI-99.0 

M-1 

IF(L.NE.l)  M-IC0UNT(ISECT,L-1)+1 

DO  10  I-M,ICOUNT( ISECT, L) 

READ(KK,END-20,ERR=20)  A,B 

IF(A.LT.XA(1) .0R.A.GT.XA(2))  GO  TO  10 

IF(B.LT.YA(1) .0R.B.GT.YA(2))  GO  TO  10 

XMA-AMAX1(A,XMA) 

XMI-AMIN1(A,XMI) 

CALL  P0INT(A,B,1) 
10   CONTINUE 
20   CONTINUE 

YMA-YA(2)+0.05*(YA(2) -YA(1)) 

CALL  SETLCH(XMA,YMA,1,0,1,I) 

WRITE(100,930)  ROTF(ISECT,L) 

PRINT  910, XMA,  YMA 

XMAX-AMAXl ( XMA , XMAX ) 
30  CONTINUE 

XC-XMAX+AMAX1( 2 . 0*RAD . 0 . 2 5* ( XA ( 2 ) - XA ( 1 ) ) +RAD ) 

YC-YA(1)+0.33*(YA(2)-YA(1)) 

CALL  CIRCLE (XC,YC, RAD) 

XAA-XC 

YAA-YC+RAD+AMAX 1 ( RAD , 0 . 1 * ( YA ( 2 ) -  YA  ( 1 ) )  ) 
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CALL  SETLCH(XAA,YAA,1,0,1,1) 

WRITE(100,920)  RAD 

CALL  FRAME 

RETURN 
900   FORMAT (//9X,8H    XMIN,8H    XEND,8H    YMIN , 8H    YEND/ 

1         9X,4F8.4//) 
910   FORMAT ( 9X, 'XMA,  YMA',2F15.5) 
920   FORMAT ( IX, 'RHO  -  ',F7.4) 
930   FORMAT (IX, 'IOTA  -  ',F7.5) 

END 


SUBROUTINE  CIRCLE (XC,YC.R) 

PLOT  THE  CIRCLE  CENTERED  AT  (XC,YC)  WITH  RADIUS  R 
PI-3. 1415926536 
PI2-PI*2.0 
DO  10  1-1,101 
UP-pTo/ioo.O+d-l) 
X-X  rk^COS(UP) 
Y-YC+R*SIN(UP) 
CALL  POINT (X.Y.l) 
10   CONTINUE 
RETURN 
END 


«^'JBROUTINE  MAXIS(R,ISGN,IFAIL) 
:      DETERMINE  THE  MAGNETIC  AXIS  IN  2D 
R-0 . 80 

DO  10  K-1,100 
CALL  PSIP(F,FP,R,ISGN) 
DELT--F/FP 

IF(ABS(DELT) .LT. 0.00001)  RETURN 
R-R+DELT 
10   CONTINUE 
PRINT  900 
IFAIL-1 
RETURN 
900   FORMAT (IX, 'THE  FIXED  POINT  DOES  NOT  CONVERGE') 
END 


SUBROUTINE  PSIP(F, FP,R, ISGN) 

2D  FLUX  FUNCTION  AND  ITS  DERIVATIVE  WITH  RESPECT  TO  R 

USE  NAME2 

USE  NAME5 

USE  NAME13 

USE  NAME17 

DIMENSION  BMF(6) ,BMFP(5) ,BMFPP(4) 

F-O.O 

FP-0 . 0 

DO  10  L-1,NT 
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ZZ-ALP*L*R 

CALL  BESSEL(BMF,BMFP,BMFPP,ZZ,L) 
X1-ALP*R+1 . 0/ ( ALP*R) 
F-F-ISGN*DEL(L,L)*X1*BMF(L) 

FP-FP-ISGN*DEL(L,L)*((ALP-1.0/(ALP*R*R))*BMF(L)+X1*ALP*L*BMFP(L)) 
10   CONTINUE 

F-F+ALP*BO*R- DELMl/R 
FP-FP+DELMl/ (R*R) +ALP*BO 
RETURN 
END 


SUBROUTINE  FMAXIS(IFL,RI) 

DETERMINE  THE  POSITION  OF  THE  MAGNETIC  AXIS 

USE  NAME13 

USE  NAME17 

DIMENSION  Y(3),  WORK(3,18) 

EXTERNAL  FCN,GCN2 

CALL  MAXIS(RAX,1,IFAIL) 

IFL-0 

DV-1 . 0 

H-20.0 

Y(1)-RAX 

DO  80  1-1,20 

FINV-DV 

Y(2)-0.0 

Y(3)-0.0 

Rl-Y(l) 

Tl-0.0 

T2-0.0 

HMAX-0.0 

NODE-3 

ITER-0 
10   T2-T2+H 

ITER-ITER+1 

IF(ITER.GT.5)  GO  TO  20 

IFAIL-1 

CALL  D02CHF (Tl , T2 , NODE , Y . TOL , 1 , HMAX , FCN , GCN2 . WORK , IFAIL) 

T1-T2 

IF(IFAIL.NE.O)  GO  TO  10 

GO  TO  30 
20   PRINT  900 

STOP  79 
30   FINV-FINV+DV 

EP-SQRT(R1**2+Y(1)**2-2.0*Y(1)->R1*COS(PI2*Y(2))) 

IF(EP.LT. 0.001)  GO  TO  90 

X-Y(l)*COS(PI2*Y(2)) 

XSUM-Rl+X 

DO  70  J-1,3 

ITER-0 
40  T2-T2+H 

ITER-ITER+1 

IF(ITER.GT.5)  GO  TO  50 

IFAIL-1 

CALL  D02CHF(T1 , T2 , NODE , Y . TOL . 1 . HMAX , FCN , GCN2 , WORK . IFAIL) 

T1-T2 
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IF(IFAIL.NE.O)  GO  TO  40 

GO  TO  60 
50   PRINT  910 

STOP  80 
60   FINV-FINV+DV 

X-Y(1)*C0S(PI2*Y(2)) 

XSUM-XSUM+X 
70   CONTINUE 

Y(l)-XSUM/5.0 
80   CONTINUE 

PRINT  920 

IFL-1 

RETURN 
;      DETERMINE  THE  POSITION  OF  THE  MAGNETIC  AXIS  IN  EACH  CROSS  SECTION 
90  Y(1)-R1 

Y(2)=-0.0 

Y(3)=0.0 

RX(1)=R1 

RY(1)=0.0 

DV-1 . 0/FLOAT(NSECT) 

FINV-DV 

Tl-0.0 

T2=0.0 

HMAX-0.0 

N0DE=3 

DO  110  K-2,NSECT 
100   T2-T2+H 

IFAIL=1 

CALL  D02CHF(T1 , T2 , NODE , Y , TOL, 1 , HMAX , FCN , GCN2 , WORK , IFAIL) 

T1=T2 

IF(IFAIL.NE.O)  GO  TO  100 

FINV=FINV+DV 

RX(K)-Y(l)*COS(Y(2)*PI2) 

RY(K)-Y(1)*SIN(Y(2)*PI2) 
110   CONTINUE 

RETURN 
900   FORMAT (//9X, 'STOP  79' ) 
910   FORMAT (//9X, 'STOP  80') 
920   F0RMAT(1 OX, 'CANNOT  FIND  THE  MAGNETIC  AXIS') 

END 


SUBROUTINE  DISLD(R1,ET) 

DRIFT  SURFACE  ISLAND  WIDTH 

USE  NAMEl 

USE  NAME3 

USE  NAME4 

USE  NAME5 

USE  NAME6 

USE  NAME9 

USE  NAMEIO 

USE  NAMEll 

USE  NAME12 

USE  NAME13 

USE  NAME14 

USE  NAME15 
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USE  NAME16 

COMMON  /FV/  FINV 

EXTERNAL  FCN,GCN2 

NOR-1 

NLEFT-NOR 

DR-DRl 

DO  70  1-1.30 

R-R1+(I-1)*DR 

ITER-0 

DV-1 . 0 

H-SDELT 

FINV-0.0 

DO  10  KK-1,N0R 

EKIN(KK)-EC 

ETA(KK)-ET 

Y(4*KK-3)-R 

Y(4*KK-2)-0.0 
10  Y(4*KK-l)-0.0 

CALL  BFIELD(Y) 

DO  20  KK-1,N0R 

Y(4*KK)-SQRT(2.0*EKIN(KK))*ETA(KK)/BVAL(KK) 
20   AMAG(KK)-EKIN(KK)*(1.0-ETA(KK)*ETA(KK))/BVAL(KK) 

Tl-0.0 

T2-0.0 

HMAX-0.0 

NODE-4 

TOL-10**(-ACC)*RADL 

Y31-Y(3) 

YVl-1.0 
30   FINV-FINV+DV 

IFAIL-1 
40  T2-T2+H 

CALL  D02CHF(T1,T2, NODE. Y,T0L.1.HMAX,FCN,GCN2. WORK. IFAIL) 

YV2-Y(3)-Y31 

IF(YV2*YV1.LT.0.0)  GO  TO  50 

Y31-Y(3) 

YV1-YV2 

T1-T2 

IF(IFAIL.EQ.O)  GO  TO  60 

IF(TOL.GT.O.O)  GO  TO  40 

PRINT  900.  TOL.  IFAIL 

RETURN 
50   PRINT  910,  Y(1).Y(2) ,Y(3) 

RETURN 
60   EP-SQRT (R**2+Y( 1 ) **2 -  2 . 0*Y( 1 ) *R*COS ( PI 2*Y ( 2 ) ) ) 

ITER-ITER+1 

X-SIN(Y(2)*PI2) 

SGN-SIGN(l.O.X) 

PRINT  920.  R,  ITER,  EP,Y(3) 

IF(EP.LT.0.06)  GO  TO  80 

IF(ITER.LT.50)  GO  TO  30 
70   CONTINUE 

IF(ITER.GT.30)  RETURN 
80   IF(EP.LT.ERR)  GO  TO  110 

KTER-O 
90   R2-R+DR 

Rl-R 

SGNl-SGN 


100 


100  EP1=EP 

KTER-KTER+1 

IF(KTER.GT.20)  RETURN 

CALL  FIXPT(ITER,R2,EP,SGN,ET) 

PRINT  927,  KTER,  Rl ,  R2 ,  R,  SGN,  EP 

R-R2 

IF(EP.LT.ERR)  GO  TO  110 

X-SGN1*SGN 

IF(X.GT.O.O)  GO  TO  90 

DR-(R-R1)*0.25 

R2-(R1+R)*0.5 

GO  TO  100 
110   PRINT  930,  ITER,  R,  EP 

RFIX-R 

DR-DR2 

Rl-RFIX+DR 

KTER-0 

PRINT  940 
120  X-ABS(R2-R1) 

IF(X.LT.ERR)  GO  TO  140 

R2=R+DR 

Rl-R 
130   KTER-KTER+1 

IF(KTER.GT.7)  GO  TO  140 

CALL  WIDTH(ITER,R2,ISGN,MTER,ET) 

PRINT  950,  KTER,  R2 ,  ISGN,  MTER 

R-R2 

IF(ISGN.GT.O)  GO  TO  120 

DR-(R-R1)*0.25 

R2=(R1+R)*0.5 

GO  TO  130 
140  WD=ABS(RFIX-R1) 

PRINT  960,  RFIX,  Rl ,  WD 

RETURN 
900   FORMAT (//9X, 'TOL.IFAIL' ,E10.2,I8) 

910   F0RMAT(7X, 'TRAPPED  PARTICLES,  Y1,Y2,Y3=  ',3F10.5) 
920   FORMAT (//9X, 'R.ITER.EP' , FIO . 7 , 110 , FIO . 7 , FIO. 7) 
927   FORMAT(3X.I4,5F15.7) 

930   F0RMAT(////9X. 'DENOMINATOR  =' ,14, '  AT  R  =', FIO. 7, '   EP=',F10.7) 
940   F0RMAT(///9X, 'KTER' ,9X, 'R' ,9X, ' ISGN' ,4X, 'NTER' ) 
950   FORMAT(//10X,I2,5X,F11.8,5X,I2,6X,I3) 
960   FORMAT (////9X, 'FIXED  POINT  AT  R   =  '  , FIO . 7 

2  //9X,'LAST  SURFACE  AT  R   =  ' , FIO . 7 

3  //9X, 'WIDTH  OF  THE  ISLAND  =  ',F10.7) 
END 


SUBROUTINE  FIXPT ( ITER , R , EP , SGN , ET) 

FIND  THE  ELLIPTIC  FIXED  POINT 

USE  NAMEl 

USE  NAME 3 

USE  NAME4 

nSE  NAME5 

USE  NAME6 

USE  NAME9 

USE  NAMEIO 
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USE  NAMEll 

USE  NAME12 

USE  NAME13 

USE  NAME14 

USE  NAME15 

COMMON  /FV/  FINV 

DIMENSION  XX(500) ,YY(500) 

EXTERNAL  FCN,GCN2 

DV-1.0*ITER 

H-SDELT*ITER 

FINV-0 . 0 

DO  20  KK-1,N0R 

EKIN(KK)-EC 

ETA(KK)-ET 

Y(4*KK-3)-R 

Y(4*KK-2)-0.0 
20   Y(4*KK-l)-0.0 

CALL  BFIELD(Y) 

DO  30  KK-l.NOR 

Y(4*KK)-SQRT(2.0*EKIN(KK))*ETA(KK)/BVAL(KK) 
30   AMAG(KK)-EKIN(KK)*(1.0-ETA(KK)*ETA(KK))/BVAL(KK) 

Tl-0.0 

T2-0.0 

HMAX-0.0 

NODE-4 

TOL-10**( -ACC)*RADL 

FINV-FINV+DV 

IFAIL-1 
40   T2-T2+H 

CALL  D02CHF (Tl , T2 , NODE , Y , TOL , 1 . HMAX . FCN , GCN2 , WORK , IFAIL) 

T1-T2 

IF(IFAIL.EQ.O)  GO  TO  50 

IF(TOL.GT.O.O)  GO  TO  40 

PRINT  900,  TOL,  IFAIL 

RETURN 
50   EP-SQRT(R**2+Y(1)**2-2.0*Y(1)->R*COS(PI2*Y(2))) 

X-SIN(Y(2)*PI2) 

SGN-SIGN(1.0,X) 

RETURN 
900   F0RMAT(//9X. 'TOL. IFAIL' ,EI0.2,I8) 

END 


SUBROUTINE  WIDTH ( ITER, R. ISGN.MTER, ET) 

FIND  THE  ISLAND  WIDTH 

USE  NAMEl 

USE  NAME3 

USE  NAME4 

USE  NAMES 

USE  NAWE6 

USE  NAME9 

USE  NAMEIO 

USE  NAMEll 

USE  NAME12 

USE  NAMEl 3 

USE  NAME14 
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USE  NAME15 

COMMON  /FV/  FINV 

DIMENSION  XX(500),YY(500) 

EXTERNAL  FCN,GCN2 

ISGN-1 

SGNl-0 . 0 

DV-1 . 0*ITER 

H-SDELT*ITER 

FINV-0 . 0 

XX(1)-R 

YY(l)-0.0 

DO  10  KK-1,N0R 

EKIN(KK)-EC 

ETA(KK)-ET 

Y(4*KK-3)-R 

Y(4*KK-2)-0.0 
10  Y(4*KK-l)-0.0 

CALL  BFIELD(Y) 

DO  20  KK-1,N0R 

Y(4*KK)-SQRT(2.0*EKIN(KK))*ETA(KK)/BVAL(KK) 
20   AMAG(KK)-EKIN(KK)*(1.0-ETA(KK)*ETA(KK))/BVAL(KK) 

Tl-0.0 

T2-0.0 

HMAX-0.0 

NODE-4 

TOL-10**( -ACC)*RADL 

Yl-0.0 

DO  50  K-2,NPRT 

FINV-FINV+DV 

IFAIL-1 
30   T2=T2+H 

CALL  D02CHF(T1 . T2 , NODE , Y , TOL , 1 , HMAX , FCN , GCN2 , WORK , IFAIL) 

T1=T2 

IF( IFAIL. EQ.O)  GO  TO  40 

IF(TOL.GT.O.O)  GO  TO  30 

PRINT  900,  TOL,  IFAIL 

RETURN 
40   MTER=K 

X=SIN(Y(2)*PI2) 

XX(K)-Y(l)*COS(Y(2)*PI2) 

YY(K)=Y(1)*SIN(Y(2)*PI2) 

IF(YY(K) .LT.Yl)  GO  TO  70 

Yl-YY(K) 

EP-ABS(YY(K)) 

IF(EP.GT.0.15)  GO  TO  60 

SGN-SIGN(1.0,X) 

IF(X.LT.O.O)  GO  TO  70 

SGNl-SGN 
50   CONTINUE 
60   ISGN--1 
70   CALL  GRAPH ( XX, YY.MTER) 

RETURN 
900   F0RMAT(//9X, 'TOL, IFAIL' ,EI0.2,I8) 

END 
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SUBROUTINE  GRAPH ( XX, YY,K) 

GRAPH  FOR  THE  ISLAND 

DIMENSION  XA(2) ,YA(2) ,XX(500) ,YY(500) 

XMIN-100000 

XMAX-- 100000 

YMIN-100000 

YMAX-- 100000 

DO  10  I-1,K 

XMAX-AMAXl ( XMAX . XX ( I ) ) 

XMIN-AMIN1(XMIN,XX(I)) 

YMAX-AMAXK YMAX , YY ( I)  ) 
10      YMIN-AMIN1(YMIN,YY(I)) 

XD-XMAX-XMIN 

YD-YMAX-YMIN 

IF(XD.GT.YD)  GO  TO  40 

XA(1)-XMIN-(YD-XD)*0.5 

XA(2)-XMAX+(YD-XD)*0.5 

YA(1)-YMIN 

YA(2)-YMAX 

GO  TO  50 
40   XA(1)-XMIN 

XA(2)-XMAX 

YA(1)-YMIN-(XD-YD)*0.5 

YA(2)-YMAX+(XD-YD)*0.5 
50   CALL  MAP(XA(1),XA(2),YA(1),YA(2),0.2.0.8,0.2,0.8) 

CALL  GAXIS(XA(1)  ,  YA(1)  .XA(2)  ,  YA(1)  ,  0  ,  0  .  0  , '"F8  .  4"  ,  3  ,XA) 

CALL  GAXIS(XA(1) . YA( 1) ,XA(1) . YA(2) ,0 ,0 , 1 , "FB . 4" , 3 , YA) 

DO  60  I-l.K 
60   CALL  POINT(XX(I) .YY(I) ,1) 

CALL  FRAME 

RETURN 

END 


SUBROUTINE  ORBIT (NS , NPHI , NPSI , NMAG , NEKIN . FTAUE) 

GUIDING  CENTER  DIFFERENTIAL  EQUATIONS 

USE  NAMEl 

USE  NAME3 

USE  NAME4 

USE  NAME5 

USE  NAME6 

USE  NAME9 

USE  NAMEIO 

USE  NAMEll 

USE  NAME12 

USE  NAME13 

USE  NAME14 

USE  NAME15 

DIMENSION  FL(15) 

EXTERNAL  FCN.GCN 

TSOL1=0.0 

CALL  INIT (NS , NPHI , NPSI , NMAG , NEKIN) 

REWIND  4 

REWIND  3 

NPAR-NS*NPHI*NPSI 

NOR"NPAR*NMAG*NEKIN 
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I0-50*IC 
DO  1  1-1.15 

1  FL(I)  — 1.0 
Xl-FLOAT(NLEFT) 
FACS-0 . 5*PI 
IF(NST.LE.O)  GO  TO  4 
REWIND  4 

DO   2  1-1,20000 

READ  (4)  (PL(I,J) ,J=1,10) 

IF(PL(I,1) .LT.0.0)  GO  TO   3 

2  CONTINUE 

3  READ  (4)  T1,T2,TSE1, NODE, NLEFT.ICOL, NOR, NTER.KTER, INI 
READ  (4)  (ETA(KK),AMAG(KK),SO(KK),PSIO(KK),PHIO(KK).SJ(KK), 

1  TJ(KK),ETAO(KK),ETAHA(KK),ETAMI(KK),PSIT(KK),PHIT(KK),BVAL(KK) 

2  Y(4*KK-3) ,Y(4*KK-2) ,Y(4*KK-1) , Y(4*KK) , EKINO(KK) ,EKIN(KK) , 

3  POT(KK) ,EKMAX(KK) ,EKMIN(KK) ,ENER(KK) ,KK=1,N0R) 
TSOL1-T2/1000000. 

GO  TO  8 

4  CONTINUE 
NLOST-0 
NLEFT-NOR 

DO  5  KK-1,N0R 

TJ(KK)=0.0 

SJ(KK)-0.0 

ETA(KK)=ETAO(KK) 

ETAMA(KK)-ETA(KK) 

ETAMI(KK)=ETA(KK) 

EKIN(KK)=EKINOCKK) 

EKMAX(KK)=EKINO(KK) 

EKMIN(KK)=EKINO(KK) 

Y(4*KK-3)=R0(KK) 

Y(4*KK-2)=U0(KK) 

5  Y(4*KK-1)=PHI0(KK) 

CALL  SFCT(Y,NLEFT,SR,UR) 
PRINT  940 
DO  6  K=1,N0R 

6  PRINT  950,SR(K),UR(K) ,PHIO(K) ,SO(K),PSIO(K) ,ETAO(K) ,EKINO(K) 
CALL  BFIELD(Y) 

DO  7  KK-1,N0R 

ENER(KK)-EKIN(KK)+POT(KK)*IEL 

Y(4*KK)-SQRT(2.0*EKIN(KK))*ETA(KK)/BVAL(KK) 

7  AMAG(KK)-EKIN(KK)*(1.0-ETA(KK)*ETA(KK))/BVAL(KK) 
NTER-IOOO 

KTER-O 

Tl-O.O 

T2-0.0 

CALL  EVAL(T2,Y,N0R) 

TSOLl-0 . 0 

Xl-FLOAT(NLEFT) 

WRITE  (3)  TEMAK,TEMrN,DENER,Xl,FAV,GAV,PAV,TEAV,ETDV,TSOLl 

8  CONTINUE 
T2=T2+SDELT 

DO  9  NN-1,NC0L 

CALL  COLL(NLEFT,FTAU,FTAUE,Y) 

IC0L=IC0L+1 

9  CONTINUE 

DO  10  KK=1,NLEFT 

IF(EKIN(KK) .LT.0.04*EC)  EKIN(KK)=0 . 04*EC 
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10  SS(KK)-AMAX1(0.05,SS(KK)) 
CALL  BFIELD(Y) 

DO  11  KK-1,NLEFT 

ENER(KK)-EKIN(KK)+POT(KK)*IEL 

Y(4*KK)-SQRT(2.0*EKIN(KK))*ETA(KK:)/BVAL(KK) 

11  AMAG(KK)-EKIN(KK)*(1.0-ETA(KK)*ETA(KK))/BVAL(KK) 
IFAII^l 

TOL-10**(-ACC)*RADL 

HMAX-0.0 

N0DE-NLEFT*4 
;     ADAMS  INTEGRATION  OF  O.D.E. 

CALL  D02CHF(T1 , T2 , NODE , Y , TOL , 1 , HMAX , FCN , GCN , WORK . IFAIL) 

T2-T1 

TSOL1-T2/1000000 

CALL  EVAL(T2,Y,N0R) 

Xl-FLOAT(NLEFT) 
;      SAVE  ORBIT  DATA 

WRITE  (3)  TEMAX,TEMIN,DENER,X1,FAV,GAV,PAV,TEAV,ETDV,TS0L1 

KTER-KTER+1 

IF(KTER.LE.IO)  GO  TO  14 

PRINT  9 10 , TSOLl , TEAV , TEMAX , TEMIN . ETDV , DENER , FAV , GAV , PAV , NLEFT 

KTER-0 
14   CALL  SECOND (TIMl) 

IF(TIMl.GE.TLIM)  GO  TO  27 

IF  ( IFAIL.  EQ.  COR.  IFAIL.  EQ.  4)  GO  TO  21 

PRINT  920, IFAIL 

GO  TO  25 
21   IF(TOL.GT.O.O)  GO  TO  25 

PRINT  930, TOL 

25  IF(NLEFT.EQ.O)  GO  TO  29 

26  IF(T2.LT.TSE)  GO  TO  8 

27  CONTINUE 

DO  28  KK-1, NLEFT 
TJ(KK)-T2/1000000 
PSIT(KK)-Y(4*KK-2) 
PHIT(KK)-Y(4*KK-1) 

28  CONTINUE 

CALL  SFCT(Y,NLEFT,SJ,UR) 

29  WRITE  (3)  (FL(I) ,1=1,10) 

WRITE  (3)  T1,T2,TSE, NODE, NLEFT, ICOL, NOR, NTER.KTER, INI 

WRITE  (3)  (ETA(KK) ,AHAG(KK) ,S0(KK).PSI0(KK) ,PHIO(KK),SJ(KK) , 

1  TJ(KK) ,ETAO(KK) .ETAMA(KK) ,ETAMI(KK) ,PSIT(KK) ,PHIT(KK) ,BVAL(KK) , 

2  Y(4*KK-3) ,Y(4*KK-2) ,Y(4*KK-1) ,Y(4*KK) ,EKINO(KK) ,EKIN(KK) , 

3  POT(KK) ,EKMAX(KK) ,EKMIN(KK) ,ENER(KK) , KK-1, NOR) 
DO  30  JC-1,N0R 

YC(4*JC-3)-Y(4*JC-3) 
YC(4*JC-2)-Y(4*JC-2) 
YC(4*JC-1)-Y(4*JC-1) 
YC(4*JC)-Y(4*JC) 

30  CONTINUE 
RETURN 

910   FORMAT(2X,5F10.3,E13.3,3F10.3,3X,I5) 

920   FORMAT(///,3X,"IFAIL=",I2) 

930   F0RMAT(///,3X, "NEGATIVE  TOL=",F8.3) 

940   FORMAT (///9X , 6X , 2HSR , 6X , 2HUR . 4X , 4HPHI0 , 6X , 2HS0 , 4X , 4HPSI0 , 4X , 

1  4HETA0,3X,5HEKIN0/) 
950   FORMAT(9X,6F8.4,E8.2) 
END 
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SUBROUTINE  EVAL(TIM,Y,NOR) 

EXPONENTIAL  DECAY  FUNCTIONALS 

USE  NAMEl 

USE  NAME 3 

USE  NAME4 

USE  NAME6 

USE  NAME9 

USE  NAMEIO 

USE  NAMEll 

USE  NAME14 

CALL  SFCT(Y,NLEFT,SS,UR) 

NL-NLEFT 

I-O 

IF(I.EQ.NL)  GO  TO  2 

I-I+l 

IF(SS(I).LT.0.98)  GO  TO  1 

14-4*1 

NL4-4*NL 

Xl-Y(I4-3) 

X2-Y(I4-2) 

X3-Y(I4-1) 

X4-Y(I4) 

X5=ETA(I) 

X6-AMAG ( I ) 

X7-S0(I) 

X8=PSI0(I) 

X9-PHI0(I) 

XIO-ETAO(I) 

Xll-ETAMA(I) 

X12=-ETAMI(I) 

XI3-EKIN0(I) 

X14-EK>LAX(I) 

X15-EKMIN(I) 

X16-EKIN(I) 

X17=ENER(I) 

X18=P0T(I) 

X19=SS(I) 

Y(I4-3)-Y(NL4-3) 

Y(I4-2)-Y(NL4-2) 

Y(I4-1)-Y(NL4-1) 

Y(I4)-Y(NL4) 

ETA(I)-ETA(NL) 

AMAG(I)-AMAG(NL) 

SO(I)-SO(NL) 

TJ(I)-TJ(NL) 

SJ(I)-SJ(NL) 

PSIT(I)-PSIT(NL) 

PHIT(I)-PHIT(NL) 

PSIO(I)-PSIO(NL) 

PHIO(I)-PHIO(NL) 

ETAO(I)-ETAO(NL) 

ETAMA(I)-ETAMA(NL) 

ETAMI(I)=ETAMI(NL) 

EKINO(I)-=EKINO(NL) 

EKMAX ( I ) -EKMAX ( NL) 

EKMIN(I)-EKMIN(NL) 

EKIN(I)-EKIN(NL) 

ENER(I)-ENER(NL) 
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POT(I)-POT(NL) 

SS(I)-SS(NL) 

Y(NL4-3)-Xl 

Y(NL4-2)-X2 

Y(NL4-1)-X3 

Y(NL4)-X4 

ETA(NL)-X5 

AMAG(NL)-X6 

PSIT(NL)-X2 

PHIT(NL)-X3 

TJ(NL)-TIM/1000000 

SJ(NL)-0. 999999 

S0(NL)-X7 

PSI0(NL)-X8 

PHI0(NL)-X9 

ETA0(NL)-X10 

ETAMA(NL)-X11 

ETAMI(NL)-X12 

EKIN0(NL)-X13 

EKMAX(NL)-X14 

EKMIN(NL)-X15 

EKIN(NL)-X16 

ENER(NL)-X17 

P0T(NL)-X18 

SS(NL)-X19 

NL-NL-1 

I-I-l 

GO  TO  1 

2  NLEFT-NL 
IF(NLEFT.EQ.O)  RETURN 
CALL  BFIELD(Y) 

DO  3  I-1,NLEFT 

3  EKIN(I)-BVAL(I)*(AMAG(I)+0. 5*BVAL(I)*Y(4*I)*Y(4*I)) 
DO  4  KK-1,NLEFT 

4  PES3(KK)-EKIN(KK)/EC 
DO   5  KK-1,NLEFT 

ETA(KK)-Y(4*KK)*BVAL(KK)/SQRT(2.0*EKIN(KK)) 
IF(ABS(ETA(KK)) .GT.l.)  ETA(KK)=ETA(KK)/(ABS(ETA(KK) )+l . OE- 14) 

5  CONTINUE 

DO  11  I-1,NLEFT 

11  PES2(I)-EKIN(I)+IEL*P0T(I) 
SAV-0 . 0 

Xl-0.0 

GAV-0 . 0 

PAV-0 . 0 

ETAV-0 . 0 

TEAV-0 . 0 

TEMAX-0 . 0 

TEMIN-1000000.0 

DO  12  I-1,NLEFT 

TEAV-TEAV+PES3(I) 

TEMAX-AMAXl ( TEMAX , PES 3 ( I) ) 

TEMIN-AMINl (TEMIN . PES3 ( I ) ) 

GAV-GAV+(1.-SO(I))*COS(FACS*(SS(I)-SO(I))/(1.-SS(I)*SO(I)))/NOR 

ETAV-ETAV+ETA(I) 

SAV-SAV+SS(I) 

12  X1-X1+(PES2 ( I )  - ENERd )  )*( PES2 ( I )  - ENER(I) ) 
FAV-2.0*(NLEFT-SAV)/ FLOAT (NOR) 
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SAV-SAV/FLOAT(NLEFT) 

ETAV-ETAV/FLOAT (NLEFT) 

TEAV-TEAV/FLOAT ( NOR ) 

DENER-SQRT(X1/NLEFT)/EC 

SDV-0 . 0 

DO  13  I-l, NLEFT 

13  SDV-SDV+(SS(I)-SAV)*(SS(I)-SAV) 
SDV-SQRT ( SDV/NLEFT) 

ETDV-0 . 0 
DO  14  I-l, NLEFT 

EKMAX(I)-AMAX1(EKMAX(I) ,PES3(I)) 
EKMIN(I)-AMIN1(EKMIN(I) .PES3(I)) 
ETAMA(I)-AMAX1(ETAMA(I) ,ETA(I)) 
ETAMI(I)-AMIN1(ETAMI(I) ,ETA(I)) 

14  ETDV-ETDV+(ETAMA(I)-ETAMI(I)) 
ETDV-ETDV/FLOAT ( NLEFT ) 
NLOST-NOR- NLEFT 
PAV-FLOAT( NLEFT) /FLOAT (NOR) 
RETURN 

END 


SUBROUTINE  COLL ( NLEFT, FTAU, FTAUE.Y) 

COLLISION  OPERATOR 

USE  NAMEIO 

USE  NAMEll 

DO  1  I-l, NLEFT 

PES4(I)-TEMP(SS(I)) 

DO  2  I-l, NLEFT 

PES5(I)-1.0+(EKIN(I)/PES4(I))**1.5 

DO  3  I-l, NLEFT 

PESl(I)-RANF(l)-0.5 

PES2(I)-FTAU*(1.0-ETA(I)*ETA(I))/PES5(I) 

DO  4  I-l, NLEFT 

PES3(I)-SIGN(1.0,PES1(I)) 

DO  5  I-l, NLEFT 

ETA(I)-PES3(I)*SQRT(PES2(I))+(1.0-FTAU/PES5(I))*ETA(I) 

DO  7  1=1, NLEFT 

PESl(I)-RANF(l)-0.5 

PES2(I)-1.5*PES4(I)/PES5(I) 

PES3(I)-PES4(I)*EKIN(I)*FTAUE/PES5(I) 

DO  9  I-l, NLEFT 
EKIN(I)-EKIN(I)+2.0*SIGN(1.0,PES1(I))*SQRT(PES3(I))-2.0*FTAUE* 

1  (EKIN(I)-PES2(I))/PES5(I) 

RETURN 

END 


SUBROUTINE  BFIELD(Y) 
EVALUATION  OF  MAGNETIC  FIELD 
USE  NAMEl 
USE  NAME2 
USE  NAME3 
USE  NAMES 
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USE  NAME9 

USE  NAME13 

USE  NAMEll 

DIMENSION  BMFPP(4) ,BMFP(5) ,BMF(6) 

DO  60  II-1,NLEFT 

R-Y(II*4-3) 

U-Y(II*4-2) 

V-Y(II*4-1) 

DO  20  1-1,3 

BB(II,I)-0.0 

DO  10  J-1,3 

DBB(II,I,J)-0.0 
10   CONTINUE 
20   CONTINUE 

DO  40  K-1,NPT 

ZZ-K*ALP*R 

CALL  BESSEL( BMF . BMFP , BMFPP , ZZ , NT) 

DO  30  L-1,NT 

UVL-(L*U-K*V)*PI2 

SN-SIN(UVL) 

CS-COS(UVL) 

DBB(II , 1 , 1)-DBB(II , 1 , 1)+K*ALP*DEL(L,K)*BMFPP(L)*SN 

DBB(II,1,2)-DBB(II,1,2)+L*DEL(L,K)*BMFP(L)*CS 

DBB(II , 1 , 3)-DBB(II ,1,3) -K*ALP*DEL(L,K)*BMFP(L)*CS 

DBB(II , 2 , 2)-DBB(II , 2 , 2) -L*L*DEL(L, K)*BMF(L)*SN/ZZ 

DBB(II,2,3)-DBB(II,2,3)+L*DEL(L,K)*BMF(L)*SN/R 

DBB(II , 3 , 3)-DBB(II , 3 , 3) -K*ALP*DEL(L,K)*BMF(L)*SN 

BB(II,1)-BB(II,1)+DEL(L,K)*BMFP(L)*SN 

BB(II,2)-BB(II,2)+L*DEL(L,K)*BMF(L)*CS/ZZ 

BB(II,3)-BB(II,3)-DEL(L,K)*BMF(L)*CS 
30   CONTINUE 
40   CONTINUE 

DBB(II,2,1)-(DBB(II,1,2)-BB(II.2))/R-DELM1/(R*R) 

DBB(II,3,1)-DBB(II,1,3) 

DBB(II,3,2)-R*DBB(II,2,3) 

BB(II,3)-BB(II,3)+BO 

UP-U*PI2 

BB(II,1)-BB(II.1)+DEL10*SIN(UP) 

BB(II,2)-BB(II.2)+DEL10*COS(UP)+DELM1/R 

DBB(II , 1 . 2)-DBB( II , 1 , 2) +DEL10*COS (UP) 

DBB(II , 2 , 2)-DBB( II , 2 , 2) -DEL10*SIN(UP) 

BVAL(II)-0.0 

DO  50  1-1,3 

BVAL(II)-BVAL(II)+BB(II,I)*BB(II,I) 
50   CONTINUE 

BVAL(II)-SQRT(BVAL(II)) 
60   CONTINUE 

IF(NST.NE.l)  CALL  ELPOT(Y) 

RETURN 

END 


FUNCTION  GCN2(T,Y) 
SPECIFIED  FUNCTION  FOR  ORBITS 
USE  NAME9 
USE  NAMEll 
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USE  NAME14 

COMMON  /FV/  FINV 

GCN2-Y(3)-FINV 

RETURN 

END 


FUNCTION  GCN(T,Y) 

SPECIFIED  FUNCTION  FOR  TRANSPORT 

USE  NAME9 

USE  NAMEll 

USE  NAMEI4 

Xl-0 . 0 

CALL  SFCT(Y,NLEFT,SR,UR) 

DO  1  I-l.NLEFT 

X1-AMAX1(X1,SR(I)) 

GCN-l.O-Xl 

RETURN 

END 


SUBROUTINE 

FCN( 

X,Y 

.YPRIME) 

RIGHT  HAND 

SIDE 

OF 

O.D. 

.E. 

USE 

NAMEl 

USE 

NAME  3 

USE 

NAME4 

USE 

NAMES 

USE 

NAME  6 

USE 

NAME  9 

USE 

NAMEIO 

USE 

NAMEll 

USE 

NAME13 

DIMENSION  GRADB(3) ,GRXBB(3) 

CALL  BFIELD(Y) 

IF(NST.EQ.l)  GO  TO  50 

DO  40   L=1,NLEFT 

R=-Y(4*L-3) 

RH0-Y(4*L) 

X1-(AMAG(L)*BVAL(L)-2.0*EKIN(L))/(BVAL(L)**3) 

DO  20  1-1,3 

GRADB(I)-0.0 

DO  10  J-1,3 

GRADB(I)-GRADB(I)+BB(L,J)*DBB(L,J,I) 
10    CONTINUE 

GRADB ( I ) -GRADB ( I ) /BVAL( L) *X1 
20   CONTINUE 

GRADB(2)-GRADB(2)/R 

GRADB ( 1 ) -GRADB ( 1 ) - 1 EL*POTR( L) / ( BVAL( L) **2 ) 

GRADB(2)-GRADB(2)-IEL*POTU(L)/(BVAL(L)**2) 

GRADB(3)-GRADB(3)-IEL*POTZ(L)/(BVAL(L)**2) 

GRXBB ( 1 ) =GRADB ( 2 ) *BB ( L , 3 ) - GRADB ( 3 ) *BB ( L ,  2 ) 

GRXBB(2)=-GRADB(3)*BB(L,1)-GRADB(1)*BB(L.3) 

GRXBB(3)-GRADB(1)*BB(L,2) -GRADB(2)*BB(L, 1) 

BDGR=0 . 0 
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DO  30  1-1.3 

BDGR-BDGR+BB ( L , I ) *GRADB ( I ) 
30   CONTINUE 

YPRIME(4*L-3)-RH0*BB(L,l)+GRXBB(l) 

YPRIME(4*L-2)-(RHO*BB(L.2)+GRXBB(2))/(R*PI2) 

YPRIME(4*L-l)-(RHO*BB(L,3)+GRXBB(3))*ALP/PI2 

YPRIME(4*L)-BDGR 
40  CONTINUE 

RETURN 
50  YPRIME(1)-BB(1.1) 

YPRIME(2)-(BB(1,2)/(Y(1)*PI2)) 

YPRIME(3)-(BB(1.3)*ALP/PI2) 

RETURN 

END 


SUBROUTINE  BESSEL( BMF , BMFP , BMFPP . Z . NT) 

EVALUATION  OF  MODIFIED  BESSEL  FUNCTIONS  AND  THEIR  DERIVATIVES 

DIMENSION  BMF(6) ,BMFP(5) ,BMFPP(4) 

REAL  S18AEF,  S18AFF 

IFAIL-0 

B0-S18AEF(Z.IFAIL) 

BMF(1)-S18AFF(Z, IFAIL) 

BMF(2)-BO-2.0*BMF(l)/Z 

DO  30  I-3,NT+2 

BMF(I)-BMF(I-2)-2.0*(I-l)*BMF(I-l)/Z 
30   CONTINUE 

BOP-BMF(l) 

BMFP(I)-(BO+BMF(2))*0.5 

DO  40  I-2,NT+1 
40   BMFP(I)-(BMF(I-1)+BMF(I+1))*0.5 

BMFPP(1)-(B0P+BMFP(2))*0.5 

DO  50  1-2, NT 
50   BMFPP(I)-(BMFP(I-1)+BMFP(I+1))*0.5 

RETURN 

END 


SUBROUTINE  INIT(NS ,NPHI ,NPSI ,NMAG .NEKIN) 

INITIALIZATION  OF  VARIABLES 

USE  NAME4 

USE  NAME6 

USE  NAME8 

USE  NAMEIO 

CALL  KEEP80  (1.3) 

CALL  FR80ID 

CALL  SURF 

HPSI-l.O/NPSI 

HPHI-l.O/NPHI 

HETA-2 . 0/NMAG 

SUM-0 . 0 

HS-1. 0/50.0 

DO  3  1-1.50 

S-(I-0.5)*HS 
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X2-2.0*(1.0-S) 

SUM-SUM+HS*X2* (1.0- TE1*S**NTE1 ) **NTE2 
T0-2.0*EC/(3.0*SUM) 
DO  5  I-1,NS 
X1-(2.*FLOAT(I)-1.)/(2.*FLOAT(NS)) 

X2-1.0-SQRT(X1) 

TS-TEMP(X2) 

DO  5  M-1,NEKIN 

CALL  ENIN(NEKIN,XX,M) 

X2-TS*XX*XX 

DO  5  J-1,NPHI 

IJ-(J-1)*NV/NPHI+1 

DO  5  K-1,NPSI 

IK-(K-1)*NU/NPSI+1 

DO  5  L-1,NMAG 

KK-I+NS*(J-1+NPHI*(K-1+NPSI*(L-1+NMAG*(M-1)))) 

EKING (KK)-X2 

S0(KK)-1.0-SQRT(X1) 

PHI0(KK)-HPHI*(J-1.0) 

PSI0(KK)-HPSI*(K-1.0) 

ETAO(KK)-  -1.0+HETA*(L-0.5) 

XX- ( SURFR( IK , I J ) - RA ( I J ) ) *S0 (KK) +RA( IJ ) 

YY-(SURFZ(IK,IJ)-ZA(IJ))*SO(KK)+ZA(IJ) 

CALL  TRANSF(XX,YY,R,U) 

RO(KK)-R 

UO(KK)-U 

AMAGO (KK)  =  1 . 0  - ETAO (KK)*ETAO (KK) 

RETURN 

END 


SUBROUTINE  ENIN(N,XX,M) 
;     ENERGY  INITIALIZATION 
USE  NAME13 
X1=SQRT(PI)/(4.0*N) 
DO  7  1=1,100 
X-0.1*I 

F--0. 5*X*EXP( -X*X)+0. 5*ERF(X) - (M-0 . 5)*X1 
IF(F.GT.0.0)GO  TO  8 

7  CONTINUE 

8  CONTINUE 
IT-0 

9  F— 0.5*X*EXP(-X*X)+0.5*ERF(X)-(M-0.5)*X1 

FP-X*X*EXP(-X*X) 
DELX--F/FP 
IT-IT+1 

IF(IT.GT.20)  GO  TO  10 
X-X+DELX 

IF(ABS(DELX) .GT. 0.0001)  GO  TO  9 
XX=-X 
RETURN 
10   PRINT  900,IT,DELX 
900   F0RMAT(//,3X, "ENERGY  INTIALIZATION  DOES  NOT  CONVERGE". 

1  3X,I2.3X,E10.2) 
STOP 
END 
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FUNCTION  ERF(X) 

ERF-0 . 0 

DX-X/50.0 

DO  1  1-1.50 

X1-(I-0.5)*DX 

ERF-ERF+EXP ( - X1*X1 ) *DX 

RETURN 

END 


FUNCTION  TEMP(S) 

TEMPERATURE  PROFILE 

USE  NAME4 

USE  NAME6 

TEMP-TO* ( 1 . 0 -TE1*S**NTE1 ) **NTE2 

RETURN 

END 


SUBROUTINE  SURF 
MAGNETIC  SURFACES 
USE  NAMES 
USE  NAME13 

COMMON  /FLUX/  RR(10, 100 , 62) , SF(IO) ,RAX(62) , SJK(62 , 62) 
DIMENSION  X(300,62) ,Y(300,62) 
1   ,ISORT(300) ,THETA(300) 

DIMENSION  SP(300) ,SD(300) ,SL(300) ,SJ(300) ,F(300) .0(300) 
1   .FP(300) .FPP(300) .FPPP(300) 
NI-9 

SLC-1.0/FLOAT(NI-1) 
DO  5  I-1,NI 
SF(I)-(I-1)*SLC 
5   CONTINUE 
REWIND  7 
DO  10  K-l.NV 
READ  (7)  RAX(K).U.V 
UP-U*PI2 

RA(K)-RAX(K)*COS(UP) 
ZA(K)-RAX(K)*SIN(UP) 
10   CONTINUE 

RA(NV+1)-RA(1) 

ZA(NV+1)-ZA(1) 

DO  15  J-1,NU+1 

DO  15  K-l.NV+1 
15   RR(l.J.K)-0.0 

DO  200  IJ-2.NI 

K-O 

NJ-1 

DO  20  J-1, 12000 

READ  (7)  R.U.V 

UP-U*PI2 

K-K+1 

X(NJ,K)-R*COS(UP) 

Y(NJ.K)-R*SIN(UP) 
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IF(K.NE.NV)  GO  TO  20 
K=-0 

NJ-NJ+1 
20   CONTINUE 

IF(K.EQ.O)  NJ-NJ-1 

DO  100  K-1,NV 

A-RAX(K) 

DO  40  J-1,NJ 

ISORT(J)-J 

B-SQRT( (X(J ,K) -RA(K) )**2+(Y(J ,K) -ZA(K))**2) 

IF(A*B.GT.O.O)  GO  TO  50 

PRINT  900.  X(J,K),  Y(J,K),  RA(K) .  ZA(K) 

STOP 

50   BDOTA-RA(K)*(RA(K)-X(J,K))+ZA(K)*(ZA(K)-Y(J,K)) 

BXA-RA(K)*(Y(J,K)-ZA(K))-ZA(K)*(X(J,K)-RA(K)) 

XX-BDOTA/A/B 

IF(ABS(XX) .GT.I.O)  GO  TO  45 

THETA(J)-AC0S(XX)/PI2 

IF(BXA.GT.O.O)  THETA(J)-1.0-THETA(J) 

GO  TO  40 
45   THETA(J)-0.0 

IF(XX.LT.O.O)  THETA(J)-0.5 
40   CONTINUE 
C      SORTING 

DO  70  L-1,NJ-1 

IFLAG-O 

DO  60  M-1,NJ-L 

MM-NJ-M+1 

XX-THETA(IS0RT(MM))-THETA(IS0RT(MM-1)) 

IF(XX.GT.O.O)  GO  TO  60 

ITEMP-ISORT(MM-l) 

IS0RT(MM-1)-IS0RT(MM) 

ISORT(MM)=ITEMP 

I FLAG- 1 
60   CONTINUE 

IF(IFLAG.EQ.O)  GO  TO  80 
70   CONTINUE 
80   CONTINUE 

IS0RT(NJ+1)=NJ+1 

X(NJ  +  I,K)=X(IS0RT(1)  ,K) 

THETA(ISORT(NJ  +  1))=-1.0+THETA(ISORT(I)) 

Y(NJ+1,K)=Y(IS0RT(1) ,K) 

HU=1.0/FLOAT(NU) 

HV-1.0/FLOAT(NV) 

MJ-0 

PHI-MJ*HU 

DO  90  J-1,NJ+1 

IF(THETA(ISORT(J)) .LT.PHI)  GO  TO  90 

MJ-MJ+1 

SURFR(MJ,K)-X(ISORT(J) ,K) 

SURFZ(MJ,K)=Y(ISORT(J) ,K) 

SP(MJ)-THETA(ISORT(J)) 

PHI-MJ*HU 
90   CONTINUE 

DF-1.0/FLOAT(NU) 

NUl-NU+1 

DO  210  I=.1,NU1 
210   SD(I)=(I-1)*DF 
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DO  35  J-1,MJ 

F(J)-SURFR(J,K) 

G(J)-SURFZ(J,K) 
35   CONTINUE 

SP(MJ+1)-SP(1)+1.0 

VM-(F(2)-F(1))/(SP(2)-SP(1)) 

VN-(F(MJ)-F(MJ-1))/(SP(MJ)-SP(MJ-1)) 

CALL  SPLIF(1,MJ,SP,F,FP,FPP,FPPP,3,0. ,3,0. ,0,0.0,IND) 

CALL  INTPL( 1 , NUl , SD , SJ , 1 , MJ , SP , F , FP , FPP , FPPP ,0) 

VM-(G(2)-G(1))/(SP(2)-SP(1)) 

VN-(G(MJ)-G(MJ-1))/(SP(MJ)-SP(MJ-1)) 

CALL  SPLIF(1,MJ,SP.G,FP, FPP, FPPP, 3,0. ,3,0. ,0,0.0,IND) 

CALL  INTPL( 1 , NUl , SD , SL , 1 , MJ , SP , G , FP , FPP , FPPP , 0 ) 

SJ(NU1)-SJ(1) 

SL(NU1)-SL(1) 

IF  (K.EQ.l)  CALL  MAP(-I. 0,1. 0,-1. 0,1. 0,0. 121, 0.706, 0.253, 0.838) 

IF  (K.EQ.l)  CALL  LINE  (-1.0,0.0,1.0,0.0,0) 

IF  (K.EQ.l)  CALL  LINE  (0.0,-1.0,0.0,1.0,0) 

IF  (K.EQ.l)  CALL  POINT(RA(K) , ZA(K) , 1) 

DO  110  J-l.NU 

IF  (K.EQ.l)  CALL  LINE  (SJ ( J ) , SL( J) , SJ ( J+1) , SL( J+1) ) 

SURFR(J,K)-SJ(J) 

SURFZ(J,K)-SL(J) 

RR(IJ,J,K)-SQRT((SJ(J)-RA(K))**2+(SL(J)-ZA(K))**2) 
110   CONTINUE 

SURFR(NU1,K)-SURFR(1,K) 

SURFZ(NU1 ,K)-SURFZ( 1 ,K) 
100   CONTINUE 

DO  115  II-1,NU1 

SURFR(II,NV+1)-SURFR(II,1) 
115  SURFZ(II,NV+1)-SURFZ(II,1) 
200   CONTINUE 

CALL  FRAME 

DO  120  I-1,NI 

DO  130  K-l.NV 
130   RR(I,NU-(-l,K)-RR(I,l,K) 

DO  140  J-1,NU 
140   RR(I,J,NV+1)-RR(I,J,1) 
120  RR(I,NU+1,NVh-1)-RR(I,1,1) 

RETURN 
900   FORMAT (IX, 'ERROR  IN  SURFACE'/IX, 'X, Y,RA, ZA' ,4F10 . 5) 

END 


SUBROUTINE  SFCT(Y,NLEFT, SR.UR) 

EVALUATION  OF  RADIAL  VARIABLE 

USE  NAME6 

USE  NAME? 

USE  NAMES 

USE  NAME13 

USE  NAMEll 

USE  NAME14 

COMMON  /FLUX/  RR( 10 , 100 , 62) , SF( 10) , RAX(62 ) , SJK(62 , 62 ) 

DO  100  L-1,NLEFT 

R-Y(4*L-3) 

U-Y(4*L-2) 
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V=Y(4*L-1) 

VX-AMOD(V,1.0) 

VX-VX+1 . 0 

VX-AMOD(VX,1.0) 

DO  10  K-1,NV 

VX-VX-HV 

IF(VX.LT.O.O)  GO  TO  20 
10   CONTINUE 
20   CONTINUE 

VX-VX+HV 

AK-VX/HV 

UP-(U-VX)*PI2 

X-R*COS(UP) 

YY-R*SIN(UP) 

B-SQRT( (RA(K) -X)**2+(ZA(K) -YY)**2) 

IF(RAX(K)*B.GT.O.O)  GO  TO  40 

SR(L)-0.0 

RETURN 
40   BDOTA-(RA(K) -X)*RA(K)+(ZA(K) -YY)*ZA(K) 

BXA-(RA(K) -X)*ZA(K) - (ZA(K) -YY)*RA(K) 

XX-BDOTA/RAX(K) /B 

IF(ABS(XX) .GT.l.O)  GO  TO  45 

THETA-AC0S(XX)/PI2 

IF(BXA.GT.O.O)  THETA=1.0-THETA 

GO  TO  46 

45  THETA-O.O 
IF(XX.LT.O.O)  THETA=0.5 

46  CONTINUE 
UR(L)-THETA 
THETAX-THETA 
DO  50  J-1,NU 
THETAX-THETAX-HU 
IF(THETAX.LT.O.O)  GO  TO  60 

50   CONTINUE 

60   AJ-(THETAX+HU)/HU 

DO  70  JJ=J,J+1 

DO  70  KK=K,K+1 
70   SJK(JJ,KK)=SSS(B,JJ,KK) 

SR(L)=AJ*(1. -AK)*SJK(J+1,K)+AJ*AK*SJK(J+1,K+1) 
1        +(1.-AJ)*(1. -AK)*SJK(J,K)+(1. -AJ)*AK*SJK(J,K+1) 
100   CONTINUE 

RETURN 

END 


FUNCTION  SSS(B,J,K) 

RADIAL  FUNCTION 

USE  NAME6 

USE  NAME7 

USE  NAMES 

USE  NAME13 

USE  NA^^E11 

USE  NAaME14 

COMMON  /FLUX/  RR( 10 , 100 , 62) , SF(IO) ,RAX(62) , SJK(62 , 62) 

DO  30  I-1,NI 

IF(B.GT.RR(I,J,K))  GO  TO  30 
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IL-I 

GO  TO  80 
30   CONTINUE 

SSS-1.0 

RETURN 
80   SSS-SF(IL-1)  +  (B-RR(IL-1,J,K))/(RR(IL,J,K)-EIR(IL-1,J,K))* 
1     (SF(IL)-SF(IL-l)) 

RETURN 

END 


SUBROUTINE  TRANSF(X, Y.R.U) 

TRANFORMATION  INTO  CYLINDRICAL  COORDINATES 

USE  NAME5 

USE  NAME13 

R-SQRT(X*X+Y*Y) 

Xl-X/R 

IF(ABS(X1) .LT.1.0)  GO  TO  5 

U-O.O 

IF(X.LT.O.O)  U-0.5 

RETURN 

U=AC0S(X1)/PI2 

IF(Y.LT.O.O)  U-l.O-U 

RETURN 

END 


SUBROUTINE  SPLIF  (M .N , S , F, FP . FPP , FPPP ,KM, VM.KN.VN .MODE, FQM, IND) 

GENERAL  PURPOSE  SPLINE  FIT 

DIMENSION  S(l)  ,  F(l),  FP(1),  FPP(l),  FPPP(l) 

IND-0 

K-IABS(N-M) 

IF  (K-1)  180,180,10 
10   K-(N-M)/K 

I-M 

J-M+K 

DS-S(J)-S(I) 

D-DS 

IF  (DS)  20,180,20 
20   DF-(F(J)-F(I))/DS 

IF  (KM-2)  30,40,50 
30  U-.5 

V-3.*(DF-VM)/DS 

GO  TO  80 
40  U-0 . 

V-VM 

GO  TO  80 
50  U--1. 

V--DS*VM 

GO  TO  80 
60   I-J 

J-J+K 

DS-S(J)-S(I) 

IF  (D*DS)  180,180,70 
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70   DF-(F(J)-F(I))/DS 
B-1./(DS+DS+U) 
U-B*DS 

V-B*(6.*DF-V) 
80   FP(I)-U 
FPP(I)-V 
U-(2.-U)*DS 
V-6 . *DF+DS*V 
IF  (J-N)  60,90,60 
90   IF  (KN-2)  100,110,120 
100  V-(6.*VN-V)/U 

GO  TO  130 
110  V-VN 

GO  TO  130 
120  V-(DS*VN+FPP(I))/(1.+FP(I)) 
130   B-V 

D-DS 
140   DS-S(J)-S(I) 

U-FPP(I)-FP(I)*V 
FPPP(I)-(V-U)/DS 
FPP(I)-U 

FP(I)-(F(J)-F(I))/DS-DS*(V+U+U)/6. 
V-U 
J-I 
I-I-K 

IF  (J-M)  140,150,140 
150   I-N-K 

FPPP(N)-FPPP(I) 
FPF(N)=B 

FP(N)-DF+D*(FPP(I)+B+B)/6. 
IND-1 

IF  (MODE)  180,180,160 
160   FPPP(J)-FQM 

V=-FPP(J) 
170   I-J 

J=J+K 

DS-S(J)-S(I) 

U=.FPP(J) 

FPPP(J)-FPPP(I)+.5*DS*(F(I)+F(J)-DS*DS*(U+V)/12.) 

v=u 

IF  (J-N)  170,180.170 
180  RETURN 
END 


SUBROUTINE  INTPL  (MI ,NI , SI , FI ,M,N, S , F, FP , FPP , FPPP , MODE) 
GENERAL  PURPOSE  INTERPOLATION  USING  TAYLOR  SERIES 
DIMENSION  SI(1),  FI(1),  S(l),  F(l),  FP(1) ,  FPP(l),  FPPP(l) 
K=IABS(N-M) 
K-(N-M)/K 
I=M 

MIN-MI 
NIN-NI 
D-S(N)-S(M) 

IF  (D*(SI(NI)-SI(MI)))  10,20,20 
10  MIN=NI 
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NIN-MI 
20   KI-IABS(NIN-MIN) 

IF  (KI)  40,40,30 
30   KI-(NIN-MIN)/KI 
40   II-MIN-Kl 

C-0. 

IF  (MODE)  60,60,50 
50   C-1. 
60   II-II+KI 

SS-SI(II) 
70   I-I+K 

IF  (I-N)  80,90.80 
80   IF  (D*(S(I)-SS))  70,70,90 
90  J-I 

I-I-K 

SS-SS-S(I) 

FPPPP-C*(FPPP(J)-FPPP(I))/(S(J)-S(I)) 

FF-FPPP(I)+.25*SS*FPPPP 

FF-FPP(I)+SS*FF/3. 

FF-FP(I)+.5*SS*FF 

FI(II)-F(I)+SS*FF 

IF  (II-NIN)  60,100,60 
100   RETURN 

END 


SUBROUTINE 

TPLOT(Y) 

LOSS 

;  RATES 

AND 

PLOT 

USE 

NAMEl 

USE 

NAME  3 

USE 

NAME4 

USE 

NAME  5 

USE 

NAME6 

USE 

NAME 10 

USE 

NAMEll 

USE 

NAMEl 2 

USE 

NAMEl 3 

USE 

NAME14 

OF  LOSS  RATES  AND  ORBITS 


DIMENSION  FL(15) 

DIMENSION  TAU1(5) ,TAU2(5) ,TAU3(5) ,DEV1(5) ,DEV2(5) ,DEV3(5) ,EX1(5) 

1  EX2(5) ,EX3(5) ,TINT(5) ,TEXTX(3) ,TEXTY(3) , SYM(4) ,DENP(100) , 

2  DENL(IOO) ,XDEN(100) 
NTOT-NOR-NLOST 
PRINT  320,NLOST,NTOT 
DO  1  1-1.15 
FL(I)--1.0 

IK2-0 

REWIND  3 

CALL  KEEP80  (1,3) 

CALL  FR80ID 

CALL  MAP(0. 0,14. 0,0. 0,12.0) 

IENDM-0 

IF(NSTU.LE.O)  GO  TO  4 

REWIND  4 

DO  2  1-1.7000 

READ  (4)  (PL(I.J),J=1,10) 
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lEND-I 

IF  (PL(I,1) .LT.O. )  GO  TO  3 

2  CONTINUE 

3  IEND-IEND-1 
IK- I END 
IK2-IK 

4  CONTINUE 
REWIND  3 

DO  6  1-1,7000 

IK2-IK2+1 

READ  (3)  (PL(IK2,J),J-1.10) 

IF(PL(IK2,1).LT.0.0)  GO  TO  7 

6  CONTINUE 

7  IK2-IK2-1 

8  IEND-IK2 
IF(NSTU.LE.O)  GO  TO  10 
REWIND  3 

DO  9  I-1,IEND 

WRITE  (3)  (PL(I,J),J=1,10) 

9  CONTINUE 

WRITE  (3)  (FL(I) ,1=1,10) 

WRITE  (3)  T1,T2,TSE, NODE, NLEFT,ICOL, NOR, NTER.KTER, INI 

WRITE  (3)  (ETA(KK) ,AMAG(KK) ,SO(KK) ,PSIO(KK) ,PHIO(KK),SJ(KK) ,TJ(KK) 

1  ,ETAO(KK) ,ETAMA(KK) ,ETAMI(KK) , PSIT(KK) , PHIT(KK) ,BVAL(KK), 

2  Y(4*KK-3) ,Y(4*KK-2) ,Y(4*KK-1) , Y(4*KK) , EKINO(KK) . EKIN(KK) , 

3  POT(KK) ,EKMAX(KK) ,EKMIN(KK) ,ENER(KK) ,KK-1,N0R) 
10   CONTINUE 

XD-1.5 
YD-8.0 
IN(5)-10 

I ENDM-MAXO ( I ENDM , I END ) 
IF(IEND.GT.7000)  IEND=7000 
TEXTX(1)-8H 
TEXTX(2)=8HTIME 
TEXTX(3)-8H 

CALL  SETLCH(4.50.11.0,1,0,1,0) 
WRITE  (100,540) 
XD-6 . 0 
YD-7.75 

TEXTY(1)=8H   LOG  (1 
TEXTY(2)-8H-S) 
TEXTY(3)=8H 
IN(l)-5 
NTER-IOOO 

NIC1-INT(1.0/(SDELT*SFREQ)) 
COLT-1.0/(SFREQ*100000) 
IL-IEND 
DO   18  I-1,IL 
PL(I,1)-PL(I,5) 
PL(I,2)-PL(I,6) 
PL(I,9)=PL(I,7) 
18   CONTINUE 

DO  20  I-1,IEND 

NTER-NTER+1 

PL(I,5)=ALOG(PL(I,5)) 

PL(I,6)-ALOG(PL(I,6)) 

PL(I,7)=ALOG(PL(I.7)) 

ICl-l+IC/4 
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IF(NTER.LT.ICl)  GO  TO  20 

NTER=0 

NL-N0R-ABS(PL(I,4)) 
20   CONTINUE 

IN2-IENDM-1 

Il-l+IENDM/10 

IN3-IENDM-I1 

DO  30  1-1,2 

CALL  SL0PE(I1,IN3,5,EX1(I) ,DEV1(I)) 

CALL  SL0PE(I1,IN3,7,EX2(I) ,DEV2(I)) 

CALL  SL0PE(I1,IN3,6,EX3(I) ,DEV3(I)) 

TINT(I)-(IN2-I1+1)*SDELT/1000000 

Il-(Il+IN2)/2 

IF(ABS(EX1(I)) .LT.l.E-10)  EX1(I)-1.E-10 

IF(ABS(EX2(I)).LT.L.E-10)  EX2(I)-1.E-10 

IF(ABS(EX3(I)) .LT.l.E-10)  EX3 (I)-l . E- 10 

TAU1(I)-1.0/EX1(I) 

TAU2(I)-1.0/EX2(I) 

TAU3(I)-1.0/EX3(I) 

IF(TAU3(I) .LT. 0.0001)  TAU3(I)=0. 0001 
30   CONTINUE 

TINTF-10.*TINT(1) 

TAUF=10.*TAU3(2) 

XMAX-0. 

DO  40  I-1,N0R 

XMAX=AMAX1 ( XMAX , ABS ( PHIT ( I ) ) /PI 2 ) 
40   CONTINUE 

XMEAN-XMAX*COLT/T I NTF 

PRINT  3 30, COLT, XMAX 

ITIMF=AINT(TIMF) 

PRINT  490,PL(IL,3) ,ITIMF 

PRINT  500,XMEAN 

CALL  PL0TB(4. ,2. , IN2 , 1 .TEXTX.TEXTY , SYM.XD , YD) 

YD-YD-3.5 

IN(l)-6 

TEXTY(1)=8H    LOG  C 

TEXTY(2)=8H0S 

CALL  PLOTBC-^.  ,2.  ,  IN2  ,  1  ,TEXTX,TEXTY,  SYM,XD,  YD) 

NIC1=INT ( 1 . 0/(SDELT*SFREQ) ) 

NIC2-NIC1/2 

PRINT  440 

IF(IEL.LT.O)  GO  TO  51 

MSTAU-AINT(TAUF) 

PRINT  510,MSTAU,MSTAU 

GO  TO  60 

51  TAUF-TAUF/(43.*43.) 
MSTAU-AINT(TAUF) 
PRINT  520.TAUF,MSTAU 
DO  52  1-1,2 

TAUl(I)-TAUl(I)/(43.*43. ) 
TAU2(I)-TAU2(I)/(43.-^43.) 
TAU3(I)=TAU3(I)/(43.*43.) 
EX1(I)-EX1(I)*(43.*43.) 
EX2(I)=EX2(I)*(43.*43.) 
EX3(I)-EX3(I)*(43.*43.) 

52  CONTINUE 

60   IF(NIC1.LE.IN2)  GO  TO  70 
NIC1-NIC2 
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NIC2-NIC1/2 

GO  TO  60 
70   CONTINUE 

PRINT  420 

PRINT  430 

PRINT  450 

TINT(1)-10.*TINT(1) 

TINT(2)-10.*TINT(2) 

TAU1(1)-10.*TAU1(1) 

TAU1(2)-10.*TAU1(2) 

EX1(1)-0.1*EX1(1) 

EX1(2)-0.1*EX1(2) 

PRINT  470,TINT(1).TAU1(I),EX1(1),DEV1(1) 

PRINT  460,TINT(2),TAU1(2),EX1(2).DEV1(2) 

PRINT  410 

PRINT  430 

PRINT  450 

DO  80  1-1,2 

TAU3(I)-10.*TAU3(I) 

EX3(I)-0.1*EX3(I) 

TAU2(I)-10.*TAU2(I) 

EX2(I)-0.1*EX2(I) 
80   CONTINUE 

PRINT  470,TINT(1),TAU3(1),EX3(1).DEV3(1) 

PRINT  460,TINT(2),TAU3(2).EX3(2),DEV3(2) 

PRINT  480 

PRINT  430 

PRINT  450 

PRINT  470,TINT(1; ,TAU2(1) ,EX2(1) ,DEV2(1) 

PRINT  460,TINT(2),TAU2(2),EX2(2),DEV2(2) 
C      HISTOGRAMS 
NINT=20 

NINT2=2*NINT+1 
NINT3=2*NINT 

CALL  DENSITY(NINT , NLEFT , NOR , SO , XDEN , DENP , DENL, - 1 ) 
DO  100  I-1,NINT2 
PL(I,4)-XDEN(I) 
100   PL(I,5)=DENP(I) 
CALL  FRAME 

CALL  SETLCH(6.0, 11. ,1,0,1,0) 
WRITE  (100,550) 
XD-1 . 5 
YD-8  . 0 
TEXTX(1)-8H 
TEXTX(2)-8H  S 
TEXTX(3)-8H 
TEXTY(1)-8H    INITI 
TEXTY(2)-8HAL 
TEXTY(3)=8H 
IN(l)-5 
IN(5)=4 
IN(2)-6 

CALL  PLOTS (4. ,2. ,NINT2 , 1 ,TEXTX,TEXTY, SYM.XD, YD) 
CALL  DENSITY(NINT , NLEFT , NOR , SJ , XDEN , DENP , DENL, 1) 
DO  103  I-1,NINT2 
PL(I,4)-XDEN(I) 
103   PL(I,5)-DENP(I) 

TEXTY(1)-8H    FINAL 
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TEXTY(2)-8H 
YD-YD-3.5 

CALL  PL0TB(4. ,2. ,NINT2 , 1 , TEXTX.TEXTY, SYM.XD , YD) 
XD-XD+7 . 2 
YD-YD+3 . 5 

CALL  DENSITY(NINT , NLEFT , NCR , ETAO , XDEN , DENP , DENL , - 1 ) 
DO  105  I-1,NINT2 
PL(I.4)-XDEN(I) 
105   PL(I,5)-DENP(I) 
TEXTX(1)-8H 
TEXTX(2)-8HETA 
TEXTY(1)-8H    INITI 
TEXTY(2)-8HAL 

CALL  PL0TB(4. ,2. , NINT2 , 1 .TEXTX.TEXTY, SYM.XD , YD) 
CALL  DENSITY(NINT , NLEFT , NOR , ETA , XDEN , DENP , DENL. 1 ) 
DO  107  I-1,NINT2 
PL(I.6)-DENL(I) 
PL(I.4)-XDEN(I) 

107  PL(I,5)-DENP(I) 
TEXTY(1)-8H   REMAIN 
TEXTY(2)-8HING 
YD-YD-3.5 

CALL  PL0TB(4. ,2. ,NINT2 , 1 .TEXTX.TEXTY, SYM.XD , YD) 

IN(l)-6 

YD-YD-3.5 

TEXTY(1)-8H    LOST 

TEXTY(2)-8H 

CALL  PL0TB(4. .2. .NINT2 , 1 .TEXTX .TEXTY . SYM.XD , YD) 

DO  108  I-l.NOR 

X1-PHIT(I)/PI2 

X2-PSIT(I)/PI2 

PES1(I)-AMOD(X1,1.0) 

PES2(I)-AMOD(X2.1.0) 

IF(PES1(I) .LT.O.O)  PES1(I)-PES1(I)+1.0 

108  IF(PES2(I) .LT.O.O)  PES2(I)=PES2(I)+1 . 0 

CALL  DENSITY(NINT. NLEFT. NOR. PHIO, XDEN, DENP, DENL. -1) 
RETURN 
C 
320   F0RMAT(//.9X. "NUMBER  LOST  -".15 . 12X. "NUMBER  REMAINING  -",16/) 
330   FORMAT (/6X, "COLLISION  TIME  -", F8 . 2 . 8X. "PERIODS  TRAVERSED  =" , 

1  F8.1/) 
410   FORMAT (///28X. "COSINE  TEST  FUNCTION"./) 
420   FORMAT(//30X."1-S  TEST  FUNCTION"./) 
430   FORMAT (6 IX. "STANDARD") 

440   F0RMAT(////.18X. "CONFINEMENT  TIME  TAU  AND  LOSS  RATE   1/TAU  "/) 
450   FORMAT (12X. "TIME  INTERVAL" . IIX. "TAU" , 6X, "LOSS  RATE" 

1  ,6X. "DEVIATION"./) 
460   FORMAT (5X. "HALF  INTERVAL" . F7 . 2 , 4X. FIO . 1 . 5X. FIO . 3 , 5X. FIO . 3) 
470   FORMAT ( 5X , " FULL  INTERVAL" . F7 . 2 . 4X . FIO . 1 , 5X . FIO . 3 . 5X . FIO . 3 ) 
480   FORMAT(///25X. "PARTICLE  LOSS  TEST  FUNCTION"./) 
490   F0RMAT(/.13X. "ENERCON  -"E13 . 3 , 12X. "CPU  TIME  ="I6. 

1  "S"/) 
500   FORMAT(/26X."MEAN  FREE  PATH  /  FIELD  PERIOD  =",F8.1,) 
510   FORMAT(//05X. "CONFINEMENT  TIME  -"14 . 3X. "EQUIVALENT  TO  ". 

1  14."  MS  AT  1  TESLA") 
520   FORMAT (//05X, "ELECTRON  CONFINEMENT  TIME  -"14 , 3X. "EQUIVALENT  TO 

1  ,14,"  MS  AT  1  TESLA") 
530   FORMAT (3X,17HEXPONENTIAL  DECAY) 


124 


540   FORMAT(36HEXPONENTIAL  DECAY  OF  EXPECTED  VALUES) 
550   FORMAT (lOHHISTOGRAMS) 
END 


SUBROUTINE  AVER(NINT , NLEFT . NOR , Y , X , XDEN , FS , FL , IND) 
COMPUTES  Y(X)  AVERAGED  OVER  OTHER  VARIABLES 
'  USE  NAME3 
DIMENSION  X(l) ,Y(1) ,XDEN(1) , FS(1) , FL(I) .NS(IOO) ,NL(100) 
SLIM-0.99 

IF(IND.LT.O.O)  SLIM-10.0 
XMAX— I.OEIO 
XMIN-l.OElO 
DO  1  1=1, NOR 
XMAX=AMAX1 ( XMAX , X ( I ) ) 

1  XMIN-AMIN1(XMIN,X(I)) 
Xl-0 . 0 

X2-0 . 0 

IF(XMAX.GT.O.O)  Xl-1.0 

IF(XMIN.LT.O.O)  X2=I.O 

XMIN-AINT(XMIN)-X2 

XMAX-AINT (XMAX) +X1 

DX=(XMAX-XMIN)/NINT 

X2-XMIN 

DO  5  J=1,NINT 

X1-X2 

FS(2*J-1)=0.0 

FL(2*J-1)=0.0 

NS(J)=0 

NL(J)-0 

X2=X1+DX 

XDEN(2*J-1)=X1 

XDEN(2*J)=X2 

DO  4  I=L,NOR 

IF(X1.LE.X(I).AND.X(I).LT.X2)  GO  TO  2 

GO  TO  4 

2  IF(SJ(I) .GT.SLIM)  GO  TO  3 
FS ( 2*J - 1 ) =FS ( 2*J - 1 ) +Y( I ) 
NS(J)=NS(J)+I 

GO  TO  4 

3  FL(2*J-1)-FL(2*J-1)+Y(I) 
NL(J)-NL(J)+1 

4  CONTINUE 

5  CONTINUE 

DO  6  J-I,NINT 

IF(NS(J).GT.O)FS(2*J-l)=FS(2*J-l)/NS(J) 
IF(NL(J).GT.O)FL(2*J-l)=FL(2*J-l)/NL(J) 
FS(2*J)-FS(2*J-1) 

6  FL(2*J)-FL(2*J-1) 
FS(2*NINT+1)=0.0 
FL(2*NINT+1)=0.0 

XDEN ( 2*NINT+1 ) =XDEN ( 2*NINT) 

RETURM 

END 
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SUBROUTINE  DENSITY(NINT , NLEFT .NOR.X.XDEN , FS , FL. IND) 

PARTICLE  DISTRIBUTION 

USE  NAME3 

DIMENSION  X(l) ,XDEN(I) ,FS(1) ,FL(I) 

SLIM-0.99 

IF(IND.LT.O.O)  SLIM-10.0 

XMAX— I.EIO 

XMIN-l.OElO 

DO  1  I-l.NOR 

XMAX-AMAXl ( XMAX , X ( I ) ) 

1  XMIN-AMIN1(XMIN,X(I)) 
Xl-0.0 

X2-0.0 

IF(XMAX.GT.O.O)  Xl-1 . 0 

IF(XMIN.LT.O.O)  X2-1.0 

XMIN-AINT(XMIN)-X2 

XMAX-A I NT ( XMAX ) +X1 

DX-(XMAX-XMIN)/NINT 

X2-XMIN 

DO  5  J-1,NINT 

X1-X2 

FS(2*J-l)-0.0 

FL(2*J-l)-0.0 

X2-X1+DX 

XDEN(2*J-1)-X1 

XDEN(2*J)-X2 

DO  4  I-1,N0R 

IF(X1.LE.X(I).AND.X(I).LT.X2)  GO  TO  2 

GO  TO  4 

2  IF(SJ(I) .GT.SLIM)  GO  TO  3 
FS(2*J-1)-FS(2*J-I)+1.0 
GO  TO  4 

3  FL(2*J-1)-FL(2*J-1)+1.0 

4  CONTINUE 

5  CONTINUE 

DO  6  J=1,NINT 
FS(2*J)-FS(2*J-1) 

6  FL(2*J)-FL(2*J-1) 

XDEN ( 2*NINT+1 ) -XDEN ( 2*NINT) 

FS(2*NINT+l)-0.0 

FL(2*NINT+l)-0.0 

RETURN 

END 


SUBROUTINE  SLOPE ( 1 1 . 12 , 1 3 , EX , DEV) 

ESTIMATION  OF  LOSS  RATE 

USE  NAME12 

NP-I2-I1+1 

All-FLOAT(NP) 

A12-.0.0 

A2 2-0.0 

Bl-0.0 

B2-0.0 

DO  1  1-11,12 

B1-B1+PL(I.I3) 
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B2-B2+PL(I , I3)*PL(I , 10) 

A12-A12+PL(I,10) 

A22-A22+PL(I . 10)*PL(I , 10) 

DET-A11*A22-A12*A12 

X1-(B1*A22-B2*A12)/DET 

X2-(B2*A11-B1*A12)/DET 

EX— X2 

DEV-0 . 0 

AV-0 . 0 

DO  2  1-11,12 

AV-AV+PL(I,I3) 

ERR-X1+X2*PL(I , 10) -PL(I , 13) 

DEV-DEV+ERR*ERR 

AV1-ABS(AV/NP) 

AV1=AMAX1(AV1 , 1 . OE- 10) 

DEV-SQRT(DEV/A11)/ABS(AV1) 

RETURN 

END 


SUBROUTINE  PLOTB  (XL, YL.N ,M,TEXTX,TEXTY, SYM.XD , YD) 
C      DRIVER  FOR  PLOT  ROUTINES 

C      PLOTS  PL(I,IN(J))  FOR  0<J<5  AS  A  FUNCTION  OF  PL(I,IN(5)). 
C     XL  AND  YL  ARE  THE  LENGTHS  IN  INCHES  OF  THE  X  AND  Y  AXES.  N  IS  THE 
C     NUMBER  OF  POINTS  PLOTTED  FOR  EACH  CURVE,  M  IS  THE  NUMBER  OF 
C      CURVES.  TEXTX  AND  TEXTY  ARE  THE  LABELS  FOR  THE  X  AND  Y  AXES. 
C      SYM(J)  IS  THE  LABEL  FOR  THE  J-TH  CURVE. 

USE  NAME12 

DIMENSION  TEXTX(3) ,  TEXTY(3),  XA(2) ,  XB(2),  SYM(4) 

XMIN-1000000. 

XMAX-- 1000000. 

YMAX-- 1000000. 

YMIN-1000000. 

DO  10  I-1,N 

XMAX-AMAXl ( XMAX , PL ( I , IN ( 5 ) ) ) 

XMIN-AMIN1(XMIN,PL(I,IN(5))) 

DO  10  J=1,M 

YMAX-AMAXl ( YMAX , PL ( I , IN ( J ) ) ) 
10   YMIN=AMIN1(YMIN,PL(I,IN(J))) 

IF  (XMAX-XMIN.LE. 0.0000001)  XMAX=XMIN+0 . 0000001 

IF  (YMAX-YMIN.LE.  0.0000001)  YMAX=\'MIN+0  .  0000001 

XA(1)-XMIN 

XA(2)-XMAX 

XB(1)-YMIN 

XB(2)-YMAX 

CALL  GAXIS  (XD,YD,XD+XL,YD,0,0,0,"F8.2" ,3,XA) 

CALL  GAXIS  (XD,YD,XD,YD+YL,0,0,1,"E8.1",3,XB) 

CALL  SETLCH  (XD+0 . 5 , YD-0 . 5 , 1 , 0 , 1 ,0) 

CALL  CRTBCD  (TEXTX, 3) 

CALL  SETLCH  (XD- 1 . 2 , YD+0 . 25 , 1 , 0 , 1 , 1) 

CALL  CRTBCD  (TEXTY, 3) 

X0R-=XA(1) 

YO-XB(l) 

DY=-(XB(2)-XB(1))/YL 

DX=-(XA(2)-XA(1))/XL 

DO  40  J=1,M 
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PLl-( PL( 1 , IN ( 5 ) ) -XOR) /DX+XD 

PL2-(PL(1,IN(J))-Y0)/DY+YD 

CALL  SETCRT  (PL1,PL2) 

DO  30  I-1,N 

PL1-(PL(I , IN(5) ) -XOR)/DX+XD 

PL2-(PL(I , IN(J) ) -YO)/DY+YD 
30   CALL  VECTOR  (PL1,PL2) 

CALL  SETLCH  (PL1+. 1 , PL2 , 1 , 0 , 1 ,0) 

CALL  CRTBCD  (SYM(J)) 
40   CONTINUE 

RETURN 

END 


SUBROUTINE  ELPOT(Y) 
C     EVALUATION  OF  THE  ELECTRIC  POTENTIAL 

C     THE  ELECTRIC  POTENTIAL  ,POT,  IS  A  FUNCTION  OF  THE  THREE 
C     VARIABLES  ,PHI,PSI,AND  S  .REPRESENTED  IN  THE  CODE  BY  Y(4*I-1), 
C     Y(4*I-2),  AND  Y(4*I-3)  RESPECTIVELY, EACH  TERM  CONTRIBUTING 
C     TO  POT  CONTAINS  THE  FACTOR  EC-0 . 5*RADL*RADL.  FOR  EXAMPLE 
C      POT(I)-EC*EO*(l-S)**XE0  IS  THE  POTENTIAL  AS  FUNCTION  OF  S. 
C     THE  DERIVATIVES  OF  POT  WITH  RESPECT  TO  S , PHI  AND  PSI  MUST 
C     ALSO  BE  COMPUTED  HERE. 

USE  NAME4 

USE  NAMES 

USE  NAME6 

USE  NAME9 

USE  NAMEIO 

USE  NAMEII 

USE  NAME13 

EOC-EO*EC 

M-17 

N-10 

DO  10  L-1,NLEFT 

POT(L)-E0C*COS((M*Y(4*L-2)-N*Y(4*L-l))*PI2) 

POTR(L)-0.0 

POTU(L)--E0C*M*SIN((M*Y(4*L-2)-N*Y(4*L-l))*PI2)/Y(4*L-3) 

POTZ(L)-E0C*N*SIN((M*Y(4*L-2)-N*Y(4*L-l))*PI2)*ALP 
10   CONTINUE 

RETURN 

END 
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